{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example notebook for reading and plotting CM1 netcdf output\n", "\n", "## Imports\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define where your output file is and read it in\n", "\n", "I've put the output for each run into its own directory, which is a good practice so you don't accidentally overwrite them. Or even better, you can run the model itself from separate directories for each run; just copy or symlink the cm1.exe executable and other needed files into each directory" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "run_name = \"bubble_10km_noshear\"\n", "\n", "data = xr.open_dataset(run_name+\"/cm1out.nc\")\n", "\n", "#data ## this will print the metadata for your file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set up some variables to see what the times and heights are" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'zh' (zh: 40)>\n",
       "array([ 0.25    ,  0.75    ,  1.25    ,  1.75    ,  2.25    ,  2.75    ,\n",
       "        3.25    ,  3.75    ,  4.25    ,  4.75    ,  5.25    ,  5.75    ,\n",
       "        6.25    ,  6.75    ,  7.25    ,  7.75    ,  8.25    ,  8.75    ,\n",
       "        9.25    ,  9.75    , 10.250001, 10.750001, 11.250001, 11.750001,\n",
       "       12.250001, 12.750001, 13.250001, 13.750001, 14.250001, 14.750001,\n",
       "       15.250001, 15.750001, 16.25    , 16.75    , 17.25    , 17.75    ,\n",
       "       18.25    , 18.75    , 19.25    , 19.75    ], dtype=float32)\n",
       "Coordinates:\n",
       "  * zh       (zh) float32 0.25 0.75 1.25 1.75 2.25 ... 18.25 18.75 19.25 19.75\n",
       "Attributes:\n",
       "    long_name:  nominal height of scalar grid points\n",
       "    units:      km\n",
       "    axis:       Z
" ], "text/plain": [ "\n", "array([ 0.25 , 0.75 , 1.25 , 1.75 , 2.25 , 2.75 ,\n", " 3.25 , 3.75 , 4.25 , 4.75 , 5.25 , 5.75 ,\n", " 6.25 , 6.75 , 7.25 , 7.75 , 8.25 , 8.75 ,\n", " 9.25 , 9.75 , 10.250001, 10.750001, 11.250001, 11.750001,\n", " 12.250001, 12.750001, 13.250001, 13.750001, 14.250001, 14.750001,\n", " 15.250001, 15.750001, 16.25 , 16.75 , 17.25 , 17.75 ,\n", " 18.25 , 18.75 , 19.25 , 19.75 ], dtype=float32)\n", "Coordinates:\n", " * zh (zh) float32 0.25 0.75 1.25 1.75 2.25 ... 18.25 18.75 19.25 19.75\n", "Attributes:\n", " long_name: nominal height of scalar grid points\n", " units: km\n", " axis: Z" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "times = data['time']\n", "heights = data['zh']\n", "heights_f = data['zf']\n", "num_times = len(times)\n", "\n", "heights " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## We might want to just focus on a subset of the domain, set that up here, and also slice the desired variables to that subset" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "## slice the fields down that we want to use\n", "xstart = -15 ## in km\n", "xend = 15\n", "ystart = -15\n", "yend = 15\n", "\n", "## use this if you want the full domain\n", "#xstart = data['xh'][0].values\n", "#xend = data['xh'][-1].values\n", "#ystart = data['yh'][0].values\n", "#yend = data['yh'][-1].values\n", "\n", "## because the scalar variables are on xh,yh,zh, we'll use those for slicing\n", "x = data['xh'].sel(xh=slice(xstart,xend))\n", "y = data['yh'].sel(yh=slice(ystart,yend))\n", "\n", "## read in a few key variables\n", "uuu = data['uinterp'].sel(xh=slice(xstart,xend),yh=slice(ystart,yend))\n", "vvv = data['vinterp'].sel(xh=slice(xstart,xend),yh=slice(ystart,yend))\n", "www = data['winterp'].sel(xh=slice(xstart,xend),yh=slice(ystart,yend))\n", "ppert = data['prspert'].sel(xh=slice(xstart,xend),yh=slice(ystart,yend))\n", "thpert = data['thpert'].sel(xh=slice(xstart,xend),yh=slice(ystart,yend))\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Map of vertical velocity and wind vectors at a specified height" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "working on 30.0 minutes, 4.75 km\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEFCAYAAADKeq1sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAyzklEQVR4nO2deZgU1fW/38PAoMAMiIgo4IIiruxL3CL6dU0QFIMiaMSNaNRgNFGDxhCXnxuJAioK7gvgiqLGBTWC4gKDICJKWAQEBJR1FGaGmTm/P6pGm7a7p7un63b1zHmfp5+uvl33ntPV937q1K1764qqYhiGYYSXetl2wDAMw0iMCbVhGEbIMaE2DMMIOSbUhmEYIceE2jAMI+SYUBuGYYSchEItIr1FZGQyBYnIYyJyaFTaEBG5PCptHxF5PlkHRWQ/EZkjIiUi0iQi/c8iMkNEXhWRpn5aUbLlBkG27UciIp1FpKe/3UpE/pli/rR/S6p5ReQeEdk50uck840XkcIk9y0QkadjpDcWkceryfuLup1g3wcTfBezPYnI0GTKjtj/zTjl3C8i7/mvbSLS3G+Di/y0WL8/6d+Wgn9DRCQ/YvvwTJYfYWeUiEwTkVkicpaf1kREJovIByJyXYw8e4jIVBH5UER+n4bN60Rk30z4nwq5EFF/C/QGPq5KEJHdgFOBo4CJwGVZ8SwgRKRG/4ufvzPQE0BV16jqPzLgWiCo6pWquo0In6tDRPYBylR1S5JmjgfeiWH7R2CDiByYZDkJUdU/pJEtaaEWkSOAmJMfVPWPqtobOBf4SFU3+F+NUtXeqjo4Dd/SYQiQ7/v0mKp+FJCdv6jqMcCxwN/8tIuB11T1KKC3iLSJynMdcAfwa+BSEdk5FYOqeruqfl1Dv1MmGUE4SEReE5GZItIedoyYROTjiH0vF5F3oyKUY0TkDf+MvltkwSLSXUT+KyLvi8hfYhlX1a2qujkquQfwnnqzdd4Ajogq9wIRuVc8vhSRp0RkvoicLSKTRGSeiBwTy57v5z0iMl1E7vfTmorIK/7Z+1kRyReRw/1jMk1EbvKz1xORsSLyiYj8zc/bQkRe8o/LUyKSJyK7i8jbvo3n/bR9/OPwHPCXCH8GiMg1/nahiEz1t4f7tqeLyGF+2qcici/wOHApMExEXpeIqxgR6elHG9NE5GoRqScib/m/e6rEiVBT8SMiT1v/d78fcSx3FpGJfp63I455kyif/ywiZ/vfHyQij0W51Bf40P9+iH+MXxPvKutcf3uaiDT09z8JeFNELhUvApsmIqf7300F+sX63REM8/+zl+Vnxvj1d6r4giB+2xCRLiJSJCJT/Dy9/XIO9T9/JiKH+T508I/BWdX4ADAMuK+afX4HRF61/tH/DwbGyyAiXcW/OvV9GeXX4+tFZLSIfCwif63OOfGi587A6yIyTERGiEgfvw5+5Nf3BSLSX0ReFK8tHuTnHeL7+aGIHFedLVXd7m82Bhb424cDb/nbU4FfRWXrAbyrquVAEXBIlP8Jf7v4VyDiXR29EflfVudvjVDVuC+8SPYDQPDE8DE/vShin4/998eAof72OD/vEGCCnzYIGAHsAzzvp70D7OJvTwZ2T+DLe0CTiLKu9LfrAx9W+QVcAoyMyLcRaALsB6wGdgI6VfkVx86v/e0PgKbAX4FL/LR/AOcBNwN9/LR6/vtSYG+8E+BnftpI4Dh/+2q8RpQP1PfT/g2c4B+XJUB+lD87A9P87d/jXT0cBjzup7UCJvvbXwP7+9tDgMv97chjPgNoHeX3zv77n4CLo//jNPwo8t/vA072tx8FjsETmj9H2X/P/48ifd4deMnfvhU4Psqf+4EjI37reH/7ZuBuf/tu4ER/++UIW02j7B9Y9Tvi1InHgN/7208DHYE+wE1+Wjfg3qjf/hpwAF7bmY7XHnoD7/jfnwD8O0Z7Otf3MfJ1p//d0cBwv5yRCfx9H2jlbzfDq48FwExgjxi/7ULgZaAg4hgd6fu+DE948/DrdHUvdmyrI/xjtQ/wlV/O/wGzfb/6Af8PaAG86dtsFHGc/hrjeFwVYetpYC1wvv/5LaCZv30RviZF7D8zYvsW/PoR5Xvc3+4fr0Pj/ZdBvepTPXNUVUVkNp7YRSMR27Mj3vcDKqLSTonKexgwWUQAdgHa+ge9OjYC+/vbzYCqS7yd8cSmV8S+S1X1BxEpBxapaomIrPLtxWOO/77SL38/YLyf9gneH3kf8DcROROYBPwH2KiqywFEZJu//8FALxG50ffvSaA58ICI7ALsAXwGLMKrDGWRjqjqNhFZJSL7A2fgXSYfAxwhIu/5u1VUHRdVXZzgd4F3Iljll10pIo2BB0VkL/+3vhArU4p+VLEfMMvf/gTvPzsQeLjKfjwnVXWtiCAiLfEaxd9j7FYSsT3Pf18FlEZs7+JHbF/5adcAI0WkPnA7sJAd63A8qurEN3h152DgdBH5tZ//m6j9W6rq/wD8tlPF3KhydkBVn8SrI7H4M95JqWs8J/3IvlxV1/jlbfK/KhaR/wIH4XUnRvJP4DRVLY5Im+e3+zV49VJFZDs14wtVrfDb33y//lW1xXZ4x/S//r67+f7fBdwVr0BVHey3o0/Eu5LfCBQCm/Dq87KoLNtFJE9VK9hROyJJ9rfP9d9j/peZJJmuj87iKWkXvIgPYCfxLtf3xjsTVtEl4n1JgrQqPgP6qdev1pWfRb06ivAaL3iXtDP87W14fVQTRWQnPy2yPy9yO1HjjN5vMd4lE3gngUXAZlUdhnfWviNGviq+Aoar10fYC3gQGAy8pV7/2qsRvsQTrmeAP+BF4Wv9Mqf5ZfYGTo6RfzteJBBNqYjsAT/1ZZ8MrFbVXwMPkfi4JOtHFbGO25d4J7pYffHRPk8ERuFdMUUfm4XsGDgk+p9PBl73P3+uqhfjnXiv9dP29f1KRHSZXwHP+r/9GOD8qP3Xikh7v+1ECmusOvhTmnjdNu9Fve70v94PeBa4EzhDRH4Tw88duj3E78oSkTy8/2BpjDznAXeKSLtYfqofNqZAvLqX6D9aineyPdavS519v/8a43hc5X9X1a21FSj268jHwIl++vFAdP94EXCsf6LuBsxP5Gc1vz1ZPakxyUTUxXiXcS3wBAa8y42PgE/Z8YzUU0QGAytU9T0RGYIn6m8CDYEBeP1JVVwHvOg32DLgNDyx/Qn/bPkcXnfFKyJym6q+JV6f8Qy8M+hPN0lUdYaI3AM8nahPLkXG++UNAtYAt+H1+/X3f89jCfLeCoyXn0ddXIPX5fOkiJwE/MjP0WA83gAe8fOiqvPEu5M/DU+cp+JdPkbyEfCEiHRnx2j0KuB5ESkDXsET3+tF5DW8SGtlBv24A3hcRG7Ai1Kmi8gs4DEROQPvvz4pls+qei7e5fiDeJF7NFPwhPbZBP5WcQRwr7/9gHg3IpvgXVaDd+k6LolyInkFOM6PUgGewr9S8LkR70SzBvgBT7waxCnrvyIyBa/7Jm5EraqdwBs9gtft9h8R6Qwcrqpj/d3OwGtnVfxZRE7BE5KJqrosRtFr8cT6CRE5J94PToEpwLMiksx/A4Cqfi8ik4BpIlIBfA78qZqIeoKI7Ip3XG/208YDT4nI+cArqrpKRFoBl6p3Q/12vON7EzBWVUtiFRw2JPWTpWG4wb8qesOPsGJ9/xBwtf7yZnP0fgNVdVKc7xrjNdiUh2pVY7OBqm73g5B3gbNVNbrLwTCSIlRCLd5IjB3G+8ZrpBmydxveXeIq3lDV24OyZySPeCOMHsG7Mfhitv1JFRE5Cu9qqhHejcxbsuySkcOESqgNwzByFRHpBtyD1w24Fq9L9ji8rruGwP+A1sCTqroolbJzYcKLYRhGLrAKOMm/ubwY755bN7ybv1/j3dBeiCfaKWFCbRiGAYjIOf5klrRGcKg3A3ir/3E7UI43Nn0T3tyFFUAHqh9h9Evfwtz1seuuLXTvvffJthuGYeQAc+bM/l5Vd0u0TzUi/Ee80UHL8Eb1vIo3/LQ0QZ5YNvbCG/HTW3+ePVkjQi3UXbt212nTQvOcI8MwQkxhocxW1e6J9hGRASQ3pLOKH/BmTE7Am3mbUDD9ceuv4M3w/Z+IxN1fVZOO3JMZR20YhlFb+BpvYlcsWrPj7OlleKL7Cl5kXZ1I5+HNMbmpalZqKmKcCBNqwzDqDKpahDc78ReIyHi8B31VifOCFGdlnok3uapARP6ONz7/GX9iW2dgnao+mo7fJtSGYRge16jqxnQzq+pEvL7paFYD3UnuOUYxsVEfhmEYQE1EuhqOwHsUR4vqdoyHRdSGYRgBoqpxV/1JFouoDcMwQo4JtWEYRsgxoTYMwwg5JtSGYRghx4TaMAwjDmeL6HUJZhe6woTaMAwjBh+J6Ba81W79VYGyhg3PMwzDiEL8VYvH4q1N97g39TzQdRETEeqIevPmzbh6aNSqVSv55JMPndgCKCqayfLly5zZM4xkWL58GbNmfeLM3ieffMjKldELuGefF6GyA96y6CfiLQw7O4tdIKEW6q+/XsKyZV87sfXccxO44oqLndgCuOGGv/DIIw84szd+/P2sW5f2DNaUKC0tZezY0U5sAXz88QxmzvzYmb177/03ZWVlTmytWfMt48bd58QWwKOPPsj111/tzN6VV17CM8885cxeMohIg5vZcU3AO/BWUk73WdU19inMjznt1KmrfvDBp05sVVRUUFFRQX5+vhN7paWlNGjQgHr13JwrlyxZTGVlJe3bH+DE3tSpb3DCCSc7sbVgwXyaNm1G69ZtnNhz+du+/PILCgoKadOmrRN7qkppaSk77bSTE3ulpaXUr1+fvLy8GpeVzGNOk+E+Ef0OGBGVfj7wO+C3GXoiXiqEWqjtedSGYSRLJoRaRAoPg80fAk2ivlsJnA4UQQNVLa+JnVSxm4mGYRg+w2HzXvxSpAHaACcAQ71ltpxG1aHuozYMw3CFiLR5E7gwwT7X4q3VJSKxtDwwLKI2DMPw+PWA/HzqJ7hP1RQ4qrSUedu3dwJmuHLMhNoIJQUUp523mIIMemLUJZIRxGyIpgm14ZSaCHCmbZigG7mCCbWRcVyIcSZI5KeJuBEmTKiNGpEropwq0b/LhNvIJibURkrUVmGujli/28TbcIUJtVEtdVWcqyPyuJhoG0FiQm3ExMQ5NUy0jUhEpB1wPdBYVQeKyGhgZ6AZcL6q/pBKeTbhxfiJAop/ehnpY8fRUNWlqho5d2Z3Vb0Y+BRI+cEtFlHXcUxMgsUi7RyjsBCaN0+8z7p1sGlTXxEZFZE6TlXHJci1WEQm4009vytVt0yo6ygm0O6pOuYm2LWCKap6bTI7ikgLoIWqni4ig4E+wEupGAuk60NEHhGRdSIyPyKtuYhMFZFF/vsuQdg24mOX5OHA/ofaj4g0E5EHgO54jw+pJyL3A2cAs1ItL6iI+jG8Z5c8EZF2HfCOqt4uItf5n5M6Ixk1wwQhvFiUXTtR1U3AJZkqL5CIWlWn461eE0k/4HF/+3HgtCBsGz9jUVvuYP+VkQiXfdS7q+q3AKr6rYi0jLWTiAwFhgK0bbuXQ/dqD9bgcxeLsI1YhG54nqqOU9Xuqtq9RYvdsu1OTmFRWe3B/ksjEpdCvVZE9gDw39c5tF2rsUZde7H/1gC3Qj0FOM/fPg94uboMX321gDVrvg3UqSomTHiCwYP7O7EFcMklQ2q8UncqjfjbNWsoLS2tkb1U+HrZMme2NmzYQHGxOzFb+vXXzmzNX7CAjocfTvnGFU7sLVv2NVu3bnViC+Ccc87g6acfc2YvVwlqeN5E4COgg4isFJELgduBE0RkEd7SY7dXV05hYTMKCgqDcPEXdOhwEL17H+/EFsBRRx1D585d086fapT1+tSpzJ47N217qbBm7VqOO/VUJ7YAht90E7fclfIcgrQZeP75bNmyxYmtVrvvzmm//S0FBQWBR9eqyjPPPMWGDesDsxFN797Hc8ABBzmzVy2FhdCqVeJXo0bO3bJVyHOMmjTUbdu2sfPOO2fQm/hMePZZBp15phNbTz/zDAcfeCBdOnVyYm/W7Nn06NbNia1EZPqGo6oi4nTN1oxS01XIRWTQyHbtnr66TZuE+w1bvJjRq1cfpaq2FJexI5mIpFyJNMDZAwY4s9X3N7+hSRN3a42GQaTBqxOZFOtcFunajgl1DpCLN5NcNvqCgro7lM2G89UNQjc8z/gZu+NvJIvVk9qNCXVIsYZnpIqd2GsvJtQhxBqbUROs/tQ+TKhDhEVERqawelS7MKEOCdawjExjJ/7ag436CAHWmIwgyfQwvlpNYSHsuWfifdasceNLBBZRZxkTacMFVs9yGxPqLGGXpYZrrL7lLibUWcAajJEtLEDITUyoHWONxAgDVg9zCxNqh1jjMMKE1cfcwYTaEdYojDBi9TI3sOF5DrDGYIQZG76XeUSkHXA90Bi4GBiDFxivVdW/plqeRdQBYyJt5AJWTzOLqi5V1Qv97WJVHaKqvwf2kjQeLWkRdYBY5TeMHCOZCS8LFwL0FZFREanjVHVcomwicjSwUNNYrcWEOiBMpI1cw7pAUmKKql6b7M6+SJ8G/CUdY9b1EQAm0kauYnU3M4hIMxF5AOguItcALwBNgLEikvJSS6GOqDdt2kR5eTn167tx0/Waca7tFRcXO10NZevWrTTKwkKgQbJ+/Xp23XVXZ/ZKS0tp2LChM3vLV6zgm1Wr6HT4iU7sVVZWUq9e7YsXVXUTcElE0p01KS/UR2jZsqWsXr3KiS1V5corL2HVqpU1KifZiOTrZcu4evjwGtlKFlVl5OjRnH3BBYHbqqysZP369fQbOJDTBw0K3N7jTz/NP269lfvHj+fDTz4J3N6ZQ4bQb+BAtm3bFrgtgOtvuol/jRnjxBbA+Mce428jRjiJrJcuXcIVV1xMmBfYDguhjqg7duzMXnvt7cSWiHDPPQ/UqAGmUrn33WcfrvzjHykpKWGnnXZK22YylJeXs3XrVt56911mz5lDty5dArNVUlJCnzPP5ONZswAoKysjPz8/MHsrVq7kpjvuAODRsWM5olevwGxt2bKF/06fjqpy7Y03MurOOwO9Ilq+YgWLlixhsKPV3AFuufFGSkpKgOD7rPfcszX33fdwYOXXJkIdUbu+JBKRtC/V04lA9mrbNnCRBmjQoAE3Xncdcz74gE+KigK11ahRIyY8/DDNd9kFgO++/z5Qez39FcELCwsZcNppgdoqmjMHVeWCc8/l7ttvD7zbatfmzXl50iS6dOoUqJ1oIutkkJG1i7pfWwi1UOcKuXID5pCDDuLSiy4K3M6+++zDpEcfpV69eqxdty5QWz26dgVg8IABNG7cOFBbnxQVMfT88xk/Zgx5eXmB2gJo0qRJ4DaM3MCEuo7h6ublCccdx20jRgQu1M2bN2f/du246LzzArUD0LVTJ8befXetvPmViFwJRGozoe6jzgWsEsfnr8OGsWr16sDtXHrRRXTt3DlwOycdf3zgNsJKnRljncyEl4Cv3GJRt0KDDGMinRgRoU3r1oHbueIPfwjchmH1PZuYUKeJVdrw0KBBg2y7YBiB4rzrQ0SWAcVABVCuqt1d+2AYRnrUmS6QkJGtPupjVTXYcVsBYtG0UZcxsXaPdX0YhmGEnGwItQJvichsERmaBfs1wqJpw7B24JpsdH0cqaqrRaQlMFVEvlLV6VVf+uI9FKBt272y4F58rHIahpENnAu1qq7239eJyGSgJzA94vtxwDiArl2729NaDCOk1Mq+6qZNbRy1iDQWkYKqbeBEYL5LH9LFomnD+CXWLtzgOqLeHZjsT2OuD0xQ1Tcc+2AYhpFTOBVqVV0KuH0UWAawqMEw4lMru0BChg3PMwzDCDkm1NVg0bRhVI+1k2AxoTYMwwg5JtQJsCjBMJLH2ktw2POo6xBFn35Kl06dnKxOkhFSfZZ1deNfM4DrleOXLF3Kfu3aObNn1BwRqQfcDDQFZqvqozUtM9QR9ZdffsHatWuc2HryyUc4++zTfvocdHTw448/cvXw4RQXu4lCZs+Zw1XDh/PXG24I3NbUd9+lZ+/eNSskSZHetGULK1avZssPP6Qu7Glw//jxTJ8xI3A7AM+++CKnDRpE0aefOrGnqjz61FP88MMPaZeRSrspKSlhypQX07YVCAUF0KpV4pc34aWviBRFvCIfh9EPaA0IsCoTboVaqJs124XCwqZObB16aCdOOOEUJ7bAW/R14aJFXHb11Xy/fn3g9iY+/zzvf/ghd993Hw8+8kigttrtsw/9+/atWSFJRsfrN21i/5NO4ujBg9kc8Elv06ZN/P2WWzj1rLNYuGhRoLZeePllBl14IfMXLODZyZMDtVXFzKIi1n33HaWlpU7sffPNcpo128WJrQCYoqrdI17jIr7rAHwMXA5ckgljoe762GOPPdl5552d2OrSpRtdunRzYgtgn7335tXnnnNia8OGDfy4dSt77rEHq7/9lsuuvpr99t2X4489NhB7+7Vrx3VXXVXzgpIQ67YtWlBeXs7Jp5xC0w4dam4zAXfccw+t99yT20eM4ID99w/MzsuvvcZFV1zBkb/6Fb26d+dXPXpQUVEReJdVrx496NWjR6A2ImnfvgPt2wf7n2WJlUCZqqqIlGeiwFALdbaobTdFmjdvzti77+b+f/+bufPm8dqbb/KvMWPYq00bDmjfPtvu1Yj8/HwO2H9/rr7iikDtlJaW0qVjR275+98DFUxVpUP79ny3dCn16+dm87QJMLwIjBGRo4H3M1FgbtYEIy1EhC6dOtGlkzc5tLw8Iyf7rHPHTTfRcrfdArXRsGFDzuzfP1Ab4P1HBx5wQOB2jOBQ1a3AhZks04S6DpMLEVsykdlxvx1IMbXvSsgwqgh/S3WMNfbsk+5lc2Q++x+N2oQJtREKMt2naaKdXayfOrOYUBuhoEpMM9m4TaCN2oIJtREqTFyNrFJYGMoVXkyoIzCRcIyDmYS/wME0c8PINKGemWjUUlavzo5IV9k2nGCBT+awiNpwS6RQZiuirrJr0bWRI5hQG+6oEshEAr0mgw/hatUqtg9VAh25bRghxoTacEMskU4kyjWJtqvEN1b5rVrtGFGbWBs5gAm1j/WnBUh1Ih1PlNMR68iujci0SLtVkXaVSJtYGyHHbiYa7glKpOPly+bNyzqOBUCZwSJqI+MUU5BcA61OpNNd4SVehFyVHhlVxyBp/w3DESbURniINSKkJt0faXRn2LTnOk7TpqGc8GJdH0ZGqRK6agUvWoAzIdLR+WrQ5WGCbYQJi6gNt8QaiZFIpFMZrlc1oiM6orYheUaOYxF1HaG4uJgNGzbUvKAEUWqsKFRVqy8vstyq15o1qY+prto/VkQe7XMSZceNqh3enFyydCkbN250YssIL6EW6srKSid2qm4clZSUOLEH3m+79a67+HTu3OrFLAP8uHUrB/XowQV//COz58xJr5AUxGnZsq/Ztm0b02fMYPLUqYmj5nhRdLR4R71WrlwJq1ezedWqXwp89P5J+F9dd0dJSQljHniA1957L6ny0qWkpISJzz3HcX360O2YY9i6bVsgdqJRVT6ZNctZuwO3bS6XCbVQb9q0iQkTnnAiZAA333kncz77zImt79ev59aRI+l57LFcdPnlfPf994Hae/WNN1j33Xc8+tRTdD/mGHodeyxPP/NM8o0yBZGrrKzkvPMGcEDXrjwxcSJnXXUVr370UU3c/wWLKys5v7SUF8rL2ffHH/m4oiKj5VdRTAFbt27lnvvuo13Hjvzpmmt44uWXA7G1YcMG/nL99bTu0IFBF17If6dPZ/Pmzbzy+uuB2Itky5YtXD18OH0HDmSDowh+/vx53HnnzU5s5TrO+6hF5GRgFJAHPKSqt8fbt3nz5gwa9Htnvt16443ObDUtLOTVZ5+lV/fuNHZwF/nAAw6ge5cudO7YkS7+67BDDqFevWDO1ZdcMoybbhrOI08+CcAZI0bwyi23cGLbthkp/4ayMt6pqOC9igr2EqG1SEbKjUVeXh6nnHgiB3XowPJvvqFRaWkgdpo3b84/hw+n329/y8zZs5k5ezZFn35KV3+NyyApLCzk37fdxshbb3UWGB16aEcOPbSjE1u5jlOhFpE84D7gBLwl1WeJyBRVXeDSjzDQsGFDjjvmGGf2jjr8cGZNm+bEVr169Tj11NN5//3/8tRTjwJQtn07/f7+d16/5hp6N2tWo/KLKip4xl+Ytxw4Ki+Pjapk5hTwSxo2bEiH9u3pULVie4D9040bN+boI47g6COO+CmtrKwsMHvRBHXirkuISGNgDF71fFdVJ9W0TNf/Sk9gsaouVdUyYBLQz7EPhgM2bdpIz56Hc/4553Bgu3YAlJSV0eeuu5jx9ddpl6uqXBshXFVxdNMAI+psk5+fn20XjNToDzynqkOB0zJRoOuuj9bANxGfVwK9IncQkaHAUIC2bfdy55mRUdq0acuQIRdzxZCBsHo1Gxcu5JMvv+SjmTO55/332a1JEw5Io9ypFRW86/dH/yYvj9vy8+mYl5dZ5406SykNq72pXEY+QF8RGRWRPE5Vx/nbbYCqO/YZuTPrWqhjhT07dIj5P3YcQNeu3d10lhmBs0tBASf37MnJe+2V9oSWSlWuKyujZ7163JGfT+/6Ng3AyBpTVPXaON+txBPr+WSo18J1TV8JO3QltgHsaTlGUnxaWcn1+fn0z8tDanFXh5HzvAiMEZF+wJRMFOhaqGcB7UVkX2AVMBAY5NgHI0fpnpdH92w7YRjVoKo/AhdkskynNxNVtRy4HHgT+BJ4VlW/cOmDkSbVTLuO9bS52vAEutrwG4zcx/lYHFX9j6oeoKr7qeqtru3Hwh7AkyQ1eUZGZN6q7XjvVY8gTdVevPzR73FIWZTtmSGGI+xujJEaCcSpgOLUT3qRD1CKfG/VypsKnooYxhPoWH4neB51FXGF2wTacIwJteGWKgGOJJPLYVUn0iayRg5i05CMjFIVhVbbjRAtmIm6RpIlVldHmsJsfdNGmLCI2ggP0c+MTifKTlfkfdLqvjFqDZs3Vz+8/8cf3fgSiQm1kXGSjkZjrRgemR5HbL/94Qf2aNIkcbmJ0qvpn7Zo2ggb1vVhBEssIY4UyupENQYvLVnCQ59/nlo+65vOCnZ1khksojaCJTJqThRBRxK5fwx6VVbS7Z57KG3ShMuOPDK+3VjEOkmYiBshx4Tap5gCu+QNilh9z9HdD5EjQaoRzsNatmTn/HwunzyZ0kaNuOq3v42/c7xuDhNpI4cwoTbcET1WOpIkxjVX0QDoceCBTJ83j6uffpqSnXZi+ODByfsQa9swQowJteGWVEZzJLj9fvjBBzN93jwArn/4YUrKyvjn8OHJP6zJRNrIIUyo6xBbtmxh4aJFfPW///HVokV0POQQzjrjDPeOxOurjrVfHH511FEwyVs4Y/+996bN/vuzftMmWuyyS9plhoElS5dSWVnJ7i1bUlBQkNNPCbQbiZnDhDoCl/3U27ZtY0uxZysvL48Wu+4amK3y8nIuvOwynpg48ae0yy6+mN+ddlpgNqulhoL5q5NOosWuu9Kja1denzqV7sceS4tDDsmQc9mjSZMmnHneeUyfMYOddtqJlrvtxu8HDuSmG27IadE2akaoh+dt2rSJ55+f5GSxzccff4gBA/oEbqeKtevW0Wr//fm/U09laQ2WpkqG+vXrc/CBB/70+abrr2fMyJHkBbQyyptvv03Xo49O+38rpqDaV+Pd92fipCnc8a8HKSgoYMbcr3b4PigqKiro2bs3rwa0MvjuLVvy9pQpDLv0UkpKSljxzTe02HXXwEVaVXliwgS6/frXrFy1KlBbVXzxxefccUe4ViEvLvbuayd6bd3q3q9QR9RNmjShf/8znUQSHTt2CUy4YrFX27Y88eCDnNm/Pw0bNgzc3gXnnstHM2fymxNPZOj55wdqq/1++3HuWWclvX+6wtqrl7cA7JdfrqSwsDBhmZm6UsrLy2PwmWdyUIcOGSkvFg0aNOCeO+6ge5cuPPzkkwz83e8Cs1WFiPD7QYPo1qULrR11Dx1yyGG0ahXurqiwIK6Whk+Hrl2767RpRU5t1uYheqtWr3bWCNMlk9Fwbfgv165bx+4tW2bbjbRw3UddWCizVTXttSVEZNCll458+qyzrk6435gxw3jhhdFHqeqMdG2lSqgj6mxQm8dTh12kYUdxTaeh17b/zkTaABNqI8TUNtE1jHQJ9c1EwzAMw4Q6JnbZZhjpY+0n81jXhxFukpkYU0UO9MEbdRcRaQL8E2gIvKWqU5LNaxG1EU5Wr05NpKvyGEZ4uRgvOK4PrEwlo0XUcajNoz9ygmSnmRuhIte7PbZsqb7a/fADAH1FZFRE8jhVHVf1QUQOA26LyroSeAl4F3gSSHqygQm1EV5SEWvr9jDcMkVVr433pap+Duww1VlEbgA2qGqZpDiLz4Q6ARZVhwAT4Jwh16NpBzwE3CEiFwHPpZLRhNowDMMBqroGOC+dvHYzsRosSjCM6rF2Eiwm1IZhGCHHhNowjBph0XTwOBNqERkhIqtEZK7/+o0r2zXFKqJhGNnE9c3Eu1V1pGObGcFGgKSHqgb+PPHt27fToEGDQG0YsbEgxg3W9WEERmVlJQ88/HDgdkaNHRu4jc8+/5ySkpLA7RjZZfPmnyfFxnv9+KN7v1wL9eUiMk9EHhGRalYhDR8WPaTGtTfeyCdFwS78oKqMHD2aLxcuDNTO9u3bOa5PH9Z9912gdnIJaw/uyKhQi8jbIjI/xqsfMBbYD+gMfAv8K04ZQ0WkSESKvv/eGkWmmTV7Ns+88ELgdkaOHs3I0aM58IADArWzZu1a1q5bx9iHHgrUTtfOnVm6bBm9jjuO+QsWBGorFzCRdktGhVpVj1fVQ2O8XlbVtapaoaqVwHigZ5wyxqlqd1Xt3qLFbpl0LyPkagVdtHgxZw0ZwrF9+vDrI48M1NYTEybw1xtuAAhcqOd89hkAj0+cyA/+QxiCoF69evzmxBNZtnw5R5xwAm9MnRqYLYCtW7fytxEjfvp9Rt3G5aiPPSI+ng7Mry7Pxo0bqaioCM6pCFasWM677ybX+DIh1t99/z13jRrFkqVLa1xWdWzatIlzLr6YZ198kWGXXsoerVoFZmvLli18GiEugQv1vHk0aNCAnXfaiUkBXyn0OflkAH744QcWLFxIWVlZYLYaNWpEj65d6Xr00Vx21VVUVlYGZquKrVu3MuU//6G4OPFN80wFK4sW/Y/Zs2elvVp9XcLlqI87RaQzoMAy4A/VZVixYhnffLOCffbZN2DX4KWXnuO55yZw3HEnBG4LIL9BA3Zt3pz92rUL3FazZs24bcQI7ho1imuGDQvUVmFhIbf/8590PPRQvly4kP32Dfa/O/3UUznp//6Pxo0bB27rhGOPpX/fvhzfuzcXDxlC/frBNp/+ffty+dChDDj9dOrVCz6matSoERUVFc6CozZt2jJhwrt069bDib1cJtSrkHfu3E3ff3+2M3tlZWXk5+cnvX+uDdfbunUrjRo1cmbPxdA817g+hpWVlU5EOlnC3PWXiVXI+/QZ+XTv3olXIX/ppWF88IHbVcjDUwNiUK+e20aeikhDuCttLFwKDFDrRBrcH0MTaQNCLtSGYYQDE+nsYo85rSE2Y9Ewag/JrPBSFya81Eos2jBqM1a/s48JdYawymzURqxehwMT6gxildqoTVh9Dg8m1IZh/AIT6XBhQp1hrIIbuY7V4fBhQh0AVtENw8gkJtQBYWJt5CJWb4NDRHqKyDMicrv/uZ3/yOcJIvKXRHltHHWA2BhrI5cwkU5pHHVfERkVkTxOVcclyqeqM0XkWuAS//NS4AIAEXkuUV4T6oAxsTZyARPplJmiqtfG+1JEDgNui0o+J86+A4F3EhkzoXaAibURZkykM4+qfg70iU4XkWZRnwcAe6vqHYnKsz5qR1hjMMKI1Ut3iMh+wK3ASSJyoR913wPsKyL/TpTXImqHWGRthAkTabeo6hJgcFRy62TyWkTtGGscRhiwephbmFBnAWskRjax+pd7mFBnibrQWLZv3+7Ezqdz5zpbPiqXKaagTtS72ogJdRapzQ3njalTmTtvnhNbz730Eu9Om+bE1oKvvgp0tfOgqK31rK5gQh3Ftm3bnNt00YiKi4t5Y2pyq6zXlC1btnDxn/7kxBbA+x9+yBMTJzqxtWnzZoZceqmTVcEB/rdoUY1tmUgnT3ExrFmT+GULB0SxYMF8vv22mmlCGaKyspJLLhnC5s2bnNg777wzuffeu3/6HHRj+teYMc7WMFy4aBFt9tzTib3y8nJ2a9GChg0bOun+2LNVK5YuW8ay5csDtwXw3fr1jLr//rTzp1qvxo4dzbnn/i5te6kyYEAfHn/8IWf2cpVQD89r3rwFTZs2c2KrXr16jB49LuUFbtOlb98zaNNmrx3Sghy+N2L4cGf9uD26deOjd97BxQr39evXZ/KECYHbqWKvtm2ZPX26s5PeEb16cXjPninnS/fE361bT1q02C2tvOnQr98ZHHZYZ2f2cpVQC3WrVq2crvrsSqQBzjjjrJjpVQ0sCMHOy8vLeJmJqI2rkLteFVxEUj6ONbk669nzV/Ts+au086fKOeec78xWLhPqro+6jPUrGulg9aZ2EuqIuq4TZHRt1C5MoGs3FlHnANYIjURY/aj9WESdI1h0bURjAl13sIg6x7DGadTmiVJGbCyizkEsuq67mEAHy+bNUN3gqJyf8CIiA0TkCxGpFJHuUd/9TUQWi8hCETkpk3brKhZZ1R3sv67bZDqing/0Bx6MTBSRg4GBwCHAnsDbInKAqtqTdDKARdi1FxNnAzIs1Kr6JcSc6NAPmKSqpcDXIrIY6Al8lEn7dR0T7NqDCbQRiaubia2BbyI+ryTOygYiMlREikSk6Pvvv3PiXG3DLpNzF/vvjFikHFGLyNtAqxhfXa+qL8fLFiMt5oMg/CXXxwF07do9+IdF1GIsws4dTJyNRKQs1Kp6fBp2VgJtIz63Adw8Fs8wwQ4xJtBGMrjq+pgCDBSRhiKyL9AemOnItuFjl9Xhwf6LuoeI9BSRZ0Tkdv9zPRF5SkTGi8hjIhJXjzM9PO90EVkJHA68JiJvAqjqF8CzwALgDeAyG/GRPapEwoTCLXbc6zaqOhO4NiKpEVCqqhcDPwKN4+XN9KiPycDkON/dCtyaSXtGzbFukeAxYc4dtmyBkpLE+5SWAtBXREZFJI/z768BICKHAbdFZT0n6vNWb1d5AdigqnEboc1MNIAdxcREu+aYONd6pqjqtfG+VNXPgT7R6SLSLOJjV2CRqt4mIteLSGdVnRurPBNq4xeYaKeHibORCBHZD7gJOFhEFgETgatE5H6gBXBPvLwm1EZCTLTjY8JspIKqLgEGRyUPSiavCXUWUdWcWq4qljCFQbxLS0vJz88P/FhmQ5grKyudL/9lhI9Q14ANGzZQVlbmxNby5ct4++03ndgC+Pbb1Tz44L3O7E2a9BTbtm3LeLmRIxmqhGztunW89OqrGbcVj2kffMDMoqKMlhn9uyJF+rXXXmbdurUZtRePoqKZvPHGa05sASxa9D/ef/89Z/beeectli9f5sxerhJqof7mm+WsXbvGia3Jk5/l5puvd2ILvIV0jzrqGCe21q9fz5VX/oHPPvs0cFvFFDD1o7lcdMWf2FzZ2MmQtIM6dKBN65hPJKiWWIKcyNfKykr+9KehzJgxPV13U+KQQw6jZcvdndgCeP75idx88w3O7N166428+OIzzuzlKqIa3lnaXbp00+nTZzuzV15eTv36tbM3yHU3S7rHMpNdKUGdHCoqKpyv6O4K110tmWxzhYUyW1W7V79nbERkUH7+yKfz869OuF9p6TC2bx99lKrOSNdWqoRalVz339ZWkYbcOZa5cIOutoo04Lw/PHxtLpnA1X1wG7ajZBiGkS3mlpW9SlnZ1cR+jhxAGfA+wEJnXhHyPmrDMAxXqOoC2B+I9xBQgLHAGajq947cAkyoDcMwIni4tTcnZXuM7zYBjwB/j/tMjqAwoTYMw/BR1dXezO/xMb69DRiGqm517JYJtWEYxo7cXAgPwA4jkFYA/wUuzMp9PRNqwzCMCLyn2P0RuDMi9QbgJrL1eGYTasMwjF9waQN4HVgFzAHWAqdkTS9teJ5hGEYUqlouMgW4EW9d7qldNIuzAy2iNgzDiEm/erAE2IN4z4l2hUXUhmEYMVBVFZGmQDk8nlVfTKgNwzDioKpbsu0DWNeHYRhG6DGhNgzDCDkm1IZhGCHHhNowDCPkmFAbhmGEHBNqwzCMkGNCbRiGEXJMqI1AmDPH3VqXS5cu4fvvv3NmzzBcE2qh/uKLz1m9epUTWw8//ACnnXaSE1sA5577O+65587qd8wAGzdupGPHdnz22Rwn9kpKSnjllclObAEsWbKIzz//zImtyspKjjyyM6+//ooTe/PmzeWww/Zl/fr1TuyNHj2SwYP7O7EF0L//KYwff78ze7lKqFchF5HvgOXZ9iMJWgBOl+bJEOa3e3LV91zwe29V3S3bTgRBqIU6VxCRoposU58tzG/35Krvuep3bSHUXR+GYRiGCbVhGEboMaHODOOy7UCamN/uyVXfc9XvWoH1URuGYYQci6gNwzBCjgm1YRhGyDGhrgEiMkBEvhCRShHpHvXd30RksYgsFBF3M2lSRERGiMgqEZnrv36TbZ8SISIn+8d0sYhcl21/kkVElonI5/4xLsq2P/EQkUdEZJ2IzI9Iay4iU0Vkkf++SzZ9rIuYUNeM+UB/YHpkoogcDAwEDgFOBu4XkTz37iXN3ara2X/9J9vOxMM/hvcBpwAHA2f7xzpXONY/xmEej/wYXp2N5DrgHVVtD7zjfzYcYkJdA1T1S1VdGOOrfsAkVS1V1a+BxUBPt97VSnoCi1V1qaqWAZPwjrWRIVR1OrAhKrkfP6/u+jhwmkufDBPqoGgNfBPxeaWfFlYuF5F5/mVvmC9rc+24RqLAWyIyW0SGZtuZFNldVb8F8N9bZtmfOoetQl4NIvI20CrGV9er6svxssVIy9o4yES/ARgL3Izn383Av4AL3HmXEqE6rilypKquFpGWwFQR+cqPXg2jWkyoq0FVj08j20qgbcTnNsDqzHiUOsn+BhEZD7wasDs1IVTHNRVUdbX/vk5EJuN14+SKUK8VkT1U9VsR2QNYl22H6hrW9REMU4CBItJQRPYF2gMzs+xTTPyGV8XpeDdIw8osoL2I7Csi+Xg3bKdk2adqEZHGIlJQtQ2cSLiPczRTgPP87fOAeFeSRkBYRF0DROR0YAywG/CaiMxV1ZNU9QsReRZYAJQDl6lqRTZ9TcCdItIZrwthGfCHrHqTAFUtF5HLgTeBPOARVf0iy24lw+7AZBEBr81NUNU3sutSbERkItAbaCEiK4F/ALcDz4rIhcAKYED2PKyb2BRywzCMkGNdH4ZhGCHHhNowDCPkmFAbhmGEHBNqwzCMkGNCbRiGEXJMqA3DMEKOCbVhGEbI+f9jYOT187LscgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "height = 5. ## the height we want (in km); the code will find the nearest level to this\n", "this_height = np.round(data['zh'].sel(zh=height, method='nearest').values,2)\n", "\n", "#for tt in range(0,num_times): ## loop over all the available times\n", "for tt in range(6,7): ## just a specific time\n", "\n", " this_time = times[tt]\n", " this_time_minutes = this_time.dt.seconds.values/60.\n", " print(\"working on \"+str(this_time_minutes)+ \" minutes, \"+str(this_height)+\" km\")\n", " \n", " fig, ax = plt.subplots()\n", "\n", " ## color contoured vertical velocity\n", " cf = plt.contourf(x,y,www.sel(time=this_time,zh=height, method='nearest'),\n", " levels=np.arange(-20,22,2), cmap='bwr', extend='both')\n", " \n", " ## vectors every nth gridpoint\n", " skip = 3 ## if this is 10, plot vector every 10th gridpoint\n", "\n", " ## create grid for x and y points\n", " x_vec, y_vec = np.meshgrid(x[::skip],y[::skip])\n", "\n", " ## plot the vectors\n", " vec = plt.quiver(x_vec,y_vec,uuu.sel(time=this_time,zh=height, method='nearest')[::skip,::skip],\n", " vvv.sel(time=this_time,zh=height, method='nearest')[::skip,::skip],\n", " scale=10.0, scale_units='inches')\n", " \n", " ## and the key vector\n", " qk = ax.quiverkey(vec, 0.8, 0.875, 2.0, r'$2 \\frac{m}{s}$', labelpos='E',\n", " coordinates='figure', fontproperties={'size':'small'})\n", " \n", " cb = plt.colorbar(cf, ax=ax, shrink=0.9)\n", " cb.ax.tick_params(labelsize='x-small')\n", "\n", " ## add a title\n", " plt.title(run_name+\" vertical velocity (m/s) height=\"+str(this_height)+\" km time=\"+str(this_time_minutes)+ \" min\", fontsize='small')\n", " \n", " ## save the figure\n", " plt.savefig(run_name+\"/www_\"+str(this_height)+\"km_\"+\"{:04d}\".format(int(this_time_minutes))+\".png\", \n", " bbox_inches='tight',dpi=150, facecolor='white')\n", "\n", " ## uncomment this if you want to see the image in the notebook, comment it if you're doing a loop\n", " plt.show()\n", " \n", " plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## make a cross section through the storm - vertical velocity and pressure perturbation" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "working on 30.0 minutes\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAETCAYAAAAS6zytAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB9bElEQVR4nO2ddXgUVxeH3xtPSIDgCe5e3N2huBUrrsVaSr8CbaEUqhQKtIXiTqG4uwR3DW7BEjwBEiKQ5H5/7CYsYT272d0w7/PMs7uzd+49OzvzmzNn7j1XSClRUFBQULB/nGxtgIKCgoKCcSiCraCgoOAgKIKtoKCg4CAogq2goKDgICiCraCgoOAgKIKtoKCg4CCYJNhCiNpCiIlGll0ghCiRZF0PIcTgJOvyCCFWmWBDfiHEGSFEtBDCW2P9MCHEISHEJiFEOvW6k8bWaw1s3b4u1P+Dm8b7KiZuO9hwSa3bGn38qMtnE0L8oNGum5HbZRFC/GlCO58JIZpqWd9RCNHe2HqSi67fKIRwF0IsVL8fK4QIVB/rc/TUtVAIkcaCtnkLIdYKIQ4KIUYaKNtP4/0UIYSnpexI0s5+IUSAEOJwgtYIIYoJIQ4IIY4IIepr2aam+rtDQohSZrQ50xK2m4sjetgPgNrA0YQVQojMQHOgOrAMGGQTy6yEEMJi/5O6rh6AG4CUcoGU8oil6rckUsqHUsrv1R97oLbZCAYAS01oqjawV8v6VUA3E+oxm6T/SxI6AFs0Po+SUlYDsgohqumocg3wqQVN7AtsllJWB2oLIXLoKZso2FLKL6SUURa0Q5N6UsrawChgmHrdz0AvoBEwTss2PwEfA52AX01tUErZ3yxLLYQ5QlBUCLFZCHFcCFEQ3vUkhRBHNcoOFkLsSfAO1NQSQmxTXxkza1YshCgvhNirvkJ+pa1xKWWklPJFktUVgACpGgW0DaiapN5eQoi/hYrLQoglQogLQohOQojlQojzQoha2tpT2zlFfTWfrl6XTgixUQixTwixQgjhJoSoot4n+4QQCQeKkxDiHyHEMSHEKPW2mYQQ69T7ZYkQwlkIkVUIsUvdxir1ujzq/bAS0LovhBAbhBB+6vcDhBC9ddRfW33nsQEYDZQGtgohPld7bM3UdXyn9lb2qT2VBup6juvzqoy1I8k2ndX75ZgQorF6XUW1B7dPCDFcvQ9WCdUdgKbN+4QQHuptJggh6iQx6WPghMb/N1XdzrdCiD+FEEeFEP9Tf+8KOAOx6n20T/0/eEgpY4EoIURWPb/9khBisVDd9bXSaPMPofLiftbzv+v8X5I00wo4qKX5C0AOdX0B6n2XS/3dHqCFHrvXCiGyq98PFEL01FVWTRVgh/r9TqCyjno/Awqr7amlfvUWqruHdUKlHYeEEF3V7/cJ1R2EEEL8JVTn/06h/4IAgJTyjfptWiBQ/d5PSnldSvkSeCaEyKRhmxfwRkoZJqW8C/hqsV+vPgi11qnPm6VCiK3q48XLkL0WQUpp9ILKEzkICFSiuEC9/qRGmaPq1wVAP/X7WeptewD/qtd1BsYCeYBV6nW7AV/1+7VAVj22BADeGnV9oX7vAhxOsAuVtzVRY7swwBvID4QAHkCpBLt0tFNT/f4gkA74HzBAve57oDswHmimXuekfr0F5EZ1YTynXjcRqKt+Pxxoh8qrclGv+wNooN4vNwE3PfugKzBE/X4bqgNQW/21gf2A0LLvxgLN1PtgrUYZJ8BL/V6guqPxVP+Hg5Nhx0RUAnkecAfSA6fV5Q4B2TXa1zw2NG1OqE8AhxP2t4Y955L8f9XUZW+jEkVnjf+jDirvLL9GW0Jj+1+BOnr+gwggA+ADnNJos6HG/shr6v+SpI2zvD2mEv4vZ2Cf+vck/E8tgJ80tjuux+62wFfq97tQiV4DtQ2ayyJ1mR1AevX7PqjPbR11a+pBAKrzrQcwW71uPDBZ/X4y0FD9m8ap15UD/la/X6rFpoT9mBnVMXMPKJv0NwNLgEIan/2BdRqfD5Lk/MKAPiT8NvX/MEb9/ieghSlaau7igumckVJKIcQp9Y9KitB4f0rjNT8Ql2RdkyTblgTWCiFAddLnBB4ZYVMYUED9Pj0Qqn7vCQwFKmmUvSWljBBCxALXpZTRQohgtFxtNTijfr2vrj8/MFu97hgqQZgGjBJCfAIsR3ULGyalvAMghEi4LSwGVBJCjFHbtxjVCT9DCOEL+AHngOuoROW1HrvWA2uEECuAaCllmBBCW/1PUR1o+vIQFAEOJZSRUsYLIcoIVQzZFcgHZLGAHaA60e5IKWOAGCFEjBDCBdXJE6zRvi5blwJ/A8+A/VLK+CTfRyf5fF59zD5EtU+lECLBO2uMyvG4qfb2FgDBQogxUso43j2etREkpQwF0Pgd8PY4P4Nq35n7vwC8TvIbfwG+RiWigcAUIURpVBfAiwbqSmAjsE0IsRp4IlUe6U71oo0wVKL+HNU5cNvIdjQ5r34NBmI03vuicmxaCyFqotrn9wCklF10VSalfAJUE0JURBUKaQxo7qf0vNUCzd+QgIuW88sUfUjQhXs6vrc45gh2aaE6k8qg8gABPITqljcHkEmjbBlUB24ZVCKWR/2eJNsncA5oJ6V8oa4v6Ymoi5OovN5xqGJXh9Tro4DPgWVCiHZSymhA8+TQfK/vxExa7gaqMMwpVBeD68ALKeXnQvXQ6BQqwdZ2Il4B1kopD0DiLflQYIeUcroQ4g8NW/T+finlSyHEC+ALYKWe+qslqesNKg8tqV0dEn+kKqY6Um3bFeA4OvaRiXYAPAFyCyHcUYmXm5QyVi14flLKB+L9uH2izVLKh2phHIIqlJCUGCGEk4bIJf4PWsSxiJTystqWaeoLxSy1rftReceX1b8hh5TyfpLt8wgh0qvtc1f/DlAd37tQecCz9OwPQ/8LwH0hRBYp5WP151FSyk3qesqhuhOtIYRoAbRRr/dGtZ+1IqV8LYS4gkrolqi3aQB8m6ToXSllN1R3WA2BOUB9oLeuutF+3Cddn/ScugKskFKOV9viqn5dCmRPUs84VP9NvPo/fgG8Un/3UKhCtY+ADFLKBCcBKWWUEMJV7Rj5oLrgm2KjvrKGLuwWwRzBDgc2oxLmhKvfUuAIcJp3r2gVhRBdUP3pAUKIHqjEfTsqb6A9oPkkeyQqT80JeI0qdvfOAwv1zl6J6jZloxDiFynlDqGKKR9CdRVNvCpLKQ8JIaYAS4UQHc34vdqYra6vM/AQlcczUAjRRv17FujZ9idgttpzBZWntBtYLIRohOrAO69rYy2sVNvjp6f+pGwAVqg9YgCklOeEKg57BNU+HwysBv5D5cW90lKPWXZIKeOEEL+iOungrUh8CawSQrxG5QGu0WazlHKu2q6vpZTaPMrdqGKsh/UZLFRx9xD1x9zAXPWx9xI4rb4oeGlcIJYBNZJUcw+YjuoORfMhVzP1b98npQwSQpj0v6h/YwLrUYnkv1q2uQL4CSF2Apc01tcDNun7/ai8/PWoQlpIKfV52LOBJUIV696YcCekg6tqz/13A+1rshGoK4RIePi7BJiry8MWqvj7v0KIeFQXvYSOBt8A81Bp2xh12caAp5RyLapjbYt6m4Em2GcXiPcdDgUF+0cI0QHwl1JO1vJdFlTxRb3dD4UQhVEJ8hkd33dA5ZSvUN92l5JSzk5S5qSUsnySdQGonmdEmPSjdNvpDsyRUnY1YZuFwCB9NgghKgGdpZRJH3Iq2Cl2LdjqJ7M/aK6Tqm481mrvF1RPwxPYJqU0ueuPgnURQvRH5RU2le/3GEppW6wu2NZAqHq0jAI6SimDbGyOgpHYtWArKCgoKLzFEQfOKCgoKHyQKIKtoKCg4CCY00skxcmUMaPMkzu3rc2wT6Tk+s2bZMqYEV9f47qCxhu4TsfEvCY6OpJ06dK/s/61nh7hT548xtXVnTRp0gHw5s3b71xd3y378uUjvLy88fRUdRBycwMn4pExMdwJCSFPliyqCt68IfjFCzK5ueHuYtyh+jw6mnhnZzKkU9shBK9iY0mfMSOg+u0Jv+P582eAS6LNSe3WtP3ZsxDSpEmHj4/29BxuWgaU37lzm5w5c+PkpL/Hl5MRvVdlfDzBDx7g5OSEv5+fwfIfMqfOnHkqpcxsuKTj4RCCnSd3bk7u22drM+yW02fP0q5bN/p2707/Xr0Mlg/Hx6x2QkK0r791K4LOnfMxceJOcuUq/F45f/93Pw8ZUpFRo6ZQtmxV/P3Bh3AAAvfu5ZMhQzg5diwAMjiYTKNHs6NuXfw8jcsftPTFC9bfvMmKfv1UDWfL9tYIf//E3x4SArNm/UpY2FO6d3+bj0qX7evWTefYsa388svGxCqTkvR3FiuWiyVLVpM7dx69Nif8fl1IKenQoweZM2Vi6Zw5ZFRffBS0I9KmvWNrG6yFEhJJBZQtXZptq1fz3fjxXLh0yfAGBli7diVbt240quzDhyrvs3nz/uTKVdhgeSklN29eJn/+ou8KXEgIF2/coHiePImf7129iquTk9FiDarRKmceP06sg4cP6f/HHzwNC4OQkHfEsXjxcri5ueutL0HAmzTpxbVrp7hx45zRtmTL5s/DhzqucmoMiTXArPnzuXHrFuuXL1fE+gNHEexUQqGCBfn9xx/p2LMnUVH6k6MZEol9+3YTHPzugD5t3vXDh/DmzWuyZctD797jdZbTxMnpAWnS+JAune97tlwLDKSYhiCdffKE0mnTvjXAiKWQjw8hERG8DHrbU+3i7dtcOPJuQkJ/f6hWrQHDhv2k32B10+7uHnTt+h1BQRd4+FB3uXfbyM6DBwZ2iAEuXr7Md+PHs2zePNzd9V9cFFI/imCnIrp37kzJYsX46tuko4tN48GDEPz8/A0XBLZsmcvEiapsmobEGiBrVn927br1nndNSAiju3bl21atEj+fefyYMr6+xlWsxuXhQ0qkTcu5J08S6y6eJQsXb99O/Kx5kRg1qifp0kUarDckBFq1GkiDBl2Ij4/XKdqazJq1mJYt2xpte1KioqLo1KsXv40bR+GCBc2uRyH1oAh2KkIIwYwpU9iycyfrNhkalaybhw9DyJbtraLq8q5jY9/w77+/0bx5P6M19ezZo1y/rhpNntTTX/Dff8TGxSV+PnP3rkqwTaSMry9nb9xINLxY9uzcuH79nR/iQzj+/nDhwglu375mVL0hIbB37womTuwLoFW0NfdDcPA9jh8/+n4hDRv08b/vvqNY4cL0/NSSaa0VHBlFsFMZ6dKl4985c+j/+efcD9ad7kGfWGzdup9Spcro/D5BqHbt+hc/v7xkyFBVa7mkIp4tG6xZM5+goCPvFgoJIeb1az6bPx+Xx48TNzwbFkaZN29UDRq7AGWE4ExYWGITg0qU4A9N0dMwrFixIty8edmg7QmUL9+AgwfX8fDhnXf2hTaOHz/KnDnTdRfQw4YtW9i8YwczpkxBT9ZChQ8MRbBTIVUqVWLogAF07dePOA2P1RhiY2PZvHk9LupudEmFS1OgcucuSrt22vP76O5RcpnChYuqLhgaha6dPUue9OkTu++FBgXxLCaG/B4eJtnPw4eU9vZWCbb6YiCBv7ZvVxmfxMsuVKgod+5c19nzIynh4b40bdqX5csn6CyT0ISfn+GHjtoIDgmh39ChLJ0zh/Tp05u8vULqRRHsVMrIL79ECMEPv/xi0nYvXrxg+HDjkphFRUVw8+YJ4uONy4L74ME9bty4QOUS+d6uVKvbhfv3KZEtW+LnSy9fUsLLC6dHj96WM2YBSqZJw6UXL4hXp11wcXJi7MqVPHrx4r12R4wYzY8/jk5IXv9ObF3bRSc29jU1a35JxoxvC+rysjNkyEhoqLYMnrp58+YNnXr1Ykj//lStVMnwBgofFIpgp1KcnZ1ZNm8e85cuZePWrUZvlzZtWiIiIoiPj9frXQPExflx6NBSxo6typ0773Z3U3d7TiRTplj8/HKyd+9xMmfK9G5FDx9Sr3hxfmyins8iJITgqChyJPSKMOGhI4DXkyd4OzvzLEaVI188eEBpf3/Oajx4TCCjWwwXLwbSvXtdoqLeZpBNar+UkiNH/uPLLwsSG/uaevXefbCbdN+EhEDu3Hn55Zf3kgnqZfg335AubVpGDR9u0nYKHwaKYKdismbJwqpFi+g9aBDXrl83ahtXV1fSpEnDixf6k+CFhED27EUZM2Y/der0YfHiYehKJPb0aQgtW35EfPw98ucv8LYCDR7duEHeDBne1v/wIf7ahg8aQl2vv5sbIVFRiZ/LZM/OmcDAd5U14aFksRL4+eXil1/aExv75r0qw8Ie8NtvTVi37kcGD15GxowGpxsEwM3NjTJlyhsuqGbxsmVs3bmTxbNm4eSknJoK76McFamcShUq8OPo0bTu0oWICOOyfc6atZiwsHdjx7pu+52cnKhTpw/ffrub+Pg4fvqpPqdOvR104+MTxv/+14iWLbuRI0fOdx92qgVTSkndGTN4evPtBEQhr1/jr28svAH83d0J0XjoOqxwYXqU1xBPjQuGk5MT8+bNQQgn5sz5NtGzjo19zbNn9/Hw8KZ06Y/56afTFCpUVVsVwPv7aMaMRYwe/T+j7D1z7hxffvMNa//9V4lb2xmnhZBZhZDqPOs2RRHsD4C+PXpQpWJFeg4cqNML1qR69dp6PTxtEQohBM7OLrRq9Q1Llw5n8uQ2PHt2n8mTB1KzZgP69h2hvYKHDwm+fBkB+Ht7J3734PVr/Fxc3paNiDBuUeMXF0eIOiRCSAjZ0qQhKDQ0sU1NfAjH1dWVKVNW0L79MGJj33D58n5GjSrN9u1/4unpQ+PGQ3FxSZIUxQBhYU/JkCGTwXLPnj2jzaefMm3SJEoUK2ZSGwrWRQghRgC9gEHGzS9rVRTB/gAQQvD3xIncuXuX6bNVE6boyyfSuXNHDh7ckfjZmEEiCRQvXpdffz1Pzpwf8fTpHYYO/ZMRIyaSPXuSrmkaon0mOJgymTK97b4WEkJITAz+zurpDY28M9Cs29/FhQevX7/TTv2ZMwm/dUunHZ6eXhQv7sfUqYOZNq0z7duPp1On3ww19Q6a+yos7CmuroYFu/fgwbRr2ZJP2rQxWFYhZdkC8X6oJr48AwghCtnSHkWwPxA8PDyY+ttv/D1rll4vOyQE0qfPpM5kp/17Q7i5edCu3Vjq1KlG+vSZE7319/p+q9UtX8aMDCtb9p3KH0RF4W9khr53UIu7v7MzIc+fJ652cXKiRMaMnHvwQOePSQiFDBjwO0uXXqVixbbJ6gNdvHg5SpTQH8O+e+8eBw4f5sfR2uYSVrAlQgjn74EfUc2w+wvQBq7a0iZFsD8gKlesSHRMDOcCA/WW8/XNRFiYarJpU7zrBDR7WLzXvzmp4oeEUChzZhrnzfvO9yGxsfg9fYq5+IWHExIb+047ZbJk4UxCXFtLWCSB/PnT4umZ5r2eIrrQ5WW3atWNcuWq67zIhePDf2vW0KZFCyVPiB0yB2LrArnUn2sCccAhIWw2TZci2B8QQgi6tG/PvMWLtX6fICz58xelTJkqOr/XhzaBe2+dRkUxsbHkHjeOe+FvBfOh+mFjejA9HKImjxBcTUhurW6vV/HiVM6VS+tFwyT7jeDAge20aVOWR4907zQpJQv//Zcun3xiegMKVkUIkeZPVJNeavKrep2w0fBTRbA/MAb07s2S//4jTGPodlLat+9D2bLV+OWX79m7d0XiekPClfR7Te9aWzgkNCKCEZs3UyJTJnL6vI2pr7h5k1be3uaHIyIiKOPkxKv4eAJfve1bXdHPD/906Zh//DgxScRcq41mEBv7hpkzRzBqVG++/noSWbP6J20mkV27tuPs7Eyt6tWT3a6CZRkPEb2BdEnWFwGKA6swYtYJK+AQExg4EhEREZw5f55TZ84QdOcOEa9eqZaICF5FRia+T1gfGxtLtixZyO7vTw5/f7L7+5Pdz48c2bMnvmbLmjVxqHhyyZE9O40aN2fBglkMG/a254Y2QalatQXjxnXk9OndDBo0GQ8PL6Pb0Zw3IGkjT54/JzPw64YNPI+KYlb9+u9suzQ8nHEZM4KBvuD6cBKCTj4+LH38mF8Twi1A+L17LD97lm937GBIo0YM7NJFdVKGhLxjbLZs70ZN/P2Nu8MIDX3Eo0d3mTPnDEWK6J/05O+/JzF88GCL5gqJiooiOCSE4AcPCA4J4X5IyDuvwQ8e8Cw0FE8PD7y9vfFOk0a1eHuTxssr8X1aHx+KFSlCudKlKVakCK5Jpw1KxQghspYGjun4/nvgY6CdEG5SSvP7npqBQ8yaXr5sWWmPM85oivOps2c5eeYMd+/fp0TRopQrXZpCBQrg4+ODd5o0qpMh4QRRv6bx8sLZ2ZmHjx69c0LdDw5+54QLDQujdvXqfNKmDa2aNiWDxgATczgceIt27T4mMDAIN/XgFF1x2FevXjJ58mfkzVuSLl1G6hQtTWHW9Kw1Z5Q5snkzk/75hwOBgdyYOBHv0FCVWGkMK79+8yY1Tp/mvocHLkKYHRLB25vzGTPSLDiY21Wq4JQ9+1sj/f05HxvLH1u28PPgwbzJmJF4KclbsSKgii0n/E5tIxi1/fZ9+1Zz7NgWvv567jvfJY3hJ5hw/vxZ2rdvSmBgEBndYsz7jWpu3LzJynXrWLluHRcvX8bfz091sU9wADScgRz+/mTMkIHomBheaTgTCQ5Ewrqw58+5cOnSe8d0uTJlKFe6NMWLFrVbERdp056SUho/YikJA4WQdYD2esqMA3yBIVKmaGhEEWwTkFJy8vRp5i9Zwt4DB947kMuXKWMVbyQiIoLN27ezct06du7dS5UKFWjfujWtmjY1awaScHxo2bIBHTp0pXPnbjrTpyYgpSQ+Pp6rV09y61YgZcr0fscr1BUK8feH+Ph40jm9Ys7Uqfw8cyZftGpFryZN8H7x4q36aQj2uCNHePriBX8mPDBMhmDj70/J27eZXrQoNYoWfUewNZ+Krtq3j/5TplCvbl2+GjKEouXrJO0q/g6a371+HcWaNV9y8uQORo9eRtGiFbXuiwQSmu3XrxtFixZn2LARZoViNEU6OCSENi1a8Enr1tSoWtVid2MJhIeHczYwMNExOXX2LHfu3aN4kSI0a9yYHl26kCtnTou2mRySK9ilhZBn0qTRe+cTEh9P78hItiqC/T62Fuynz56xZPly5i5eTGRkJL26dqVZ48Y2uVVMrngn9L/euXMbY8Z8zeHD53jw4N1jTlfPkDt3rvDDD5+QJUsx+vSZiZdXOvz9IS4ujvDwMF6+fEbFioUJCrrG9ev7ePAghBXLF7L2338p4OmJq4uLKn1qQgNJBFsGB1Pk5EkWZcpEpbAw88UaVIIN/OLhwR1XV2bUqqVarynWGnM+hqdNy9zVq1mwaRPbdhzhwoVzHD9+FV/fTPj6ZiJDhtK4urojhEAIkWj6/v2LuHp1C8OHz8TbO2nEU0VS0ZbyPlWqfMS5czfx9fU1WrATRHrF2rU8ePiQNs2b0751a2pWq4ZzQp/1FCIiIoJTZ8+yYs0alq9eTfkyZejdrRstmza1eY+X5Aq2EOKU9PYuq6/M4/h4skZGbpNSNjG3HXNQBFsHcXFx7Ni9m3lLlrBz716aN25M727dqFmtmt3keUgq3r+OHctnffro3SZBsKWUVKnyET///AdFijR4p4y+rnwxMVH8/vtwnj9/wMSJaxk/vjN7967A2zs9mTJlZuPGCxw/HsDevcvImDETbZvUo2qlSogHD3gvzpBEsE9dvUqHy5e5njOnqnxyBBvA25s7mTNT7t49Qlq3xs3Z+X3BBpWiqt9LPz+EECzbtJvNm9dz//4Tnj9/yqRJy7lw4Q7Dh9cnbdqMpEmTibZtx9KqlWpGGX3eWFLBnjFjBDExMfz225TEdfpEOz4+nk/79GHP/v20bdGC9mpPOqVFWhdRUVGs3biRuYsXc/7CBTq3b0+vrl0pVbKkTexRBNvGpKRgB92+zdxFi1jw77/4Z8tGr65d6dSuHenSafee7IWg27ep1rAhs//8k6aNG2stk3R049KlC1iyZBlz525PXGdsv+vAwHuULJmTyMgI3N09yZ79rXgk6GCiCGkKs2YjSQT7ywMH8I6KYlymTKp1FhBs/P2p8fAhX+fMSfOSJXULtuZnjdnVNc0EuHs3hpcvn/HixVMiI30pWdK4UEBCExER4dSvn5cDB06+N5u6LtH+33ffcfzUKXasW2dz79UQQbdvM3/JEhb8+y9ZMmem16ef0uWTT1L0/EnNgm0frqKdMG3WLCrUrk3Eq1dsXb2a4wEBDOjd2+7FGiBvnjysWbKEHp99xplz78/srW0oert2nbhx4yKnTh00ub0EofLy8tYq1kajVsOnMTEsffyYzgmT7lqQrlmzMiU42Kg8Kgloiqfmb8qVy51MmfzJn/8jo8VakwUL/qBq1frviTVo/49mzpvHhi1bWLt0qd2LNaiOw3HffUdQYCA/jR5NwMGDFC5XjiPHdPW5UDAFq3XrE0LMA5oBj6WUJdTrfgeaA6+Bm0BPKeVza9lgLFJKRo8fz4q1azkREEDePHlsbZJZVK5YkX8mT6ZFx44c2bWLHNmz6y3/7Jk73333F99805N1687x4oXx3fYSSNojJAGdt/hJvWtU+3/gyZN8miULRbSkVI0HDqE6aKR6QctrZkCbW9Qza1bmPHzIjBs3+MzAPtHEh/BEEdXs1pe0y5+xXLlyjn//ncaaNaeNKr9t506+//lnDm7fnuyeQSmNs7MzjerXp1H9+mzdsYOWnToxf/p0nXd/CsZhTQ97AZD039kJlJBSfgRc4/2BRClObGwsfYcMYceePRzaudNhxTqBdq1aMaR/f5p98gnh6tGD+hI9NWjQmpIlKzJ58jcmt2WUWCcNh2hh+dWrXHjxgp80+ktrMhfogioBzy/Ab+rld2AiMAn4A6gLaMv67erkxKLChRkdGMgNPQOGDKGr66IxvHnzhlGjejB8+G9ky5ZD73B1gPMXLtC1Xz9WL1lCgfz5zbTYPmjSsCEb//uP3oMHs2DpUlub49BYTbCllPuB0CTrdkgpExI8HAWMywRvJaKiomj76afcu3+fPZs2vT8TioPyv88/p2K5cnTs2ZOwWE+tZTQF49tv/2TbtpWcP3/A6DZ0ibUpSCm5EhrK53v3srhyZTwSHuZqet/AP8AsYLd62aVedgI71Mt2oD8wU1tDDx9SxMuL0cWL023bNkIT8m7ri6mrSXqnYO5vnTnzZ7Jk8adNmx4Gy4Y8eEDzDh34c8IEqlWubF6DdkalChUI2LKFsb/8wm+TJ5sUnnJUhBBNhRCzhRDrhBD1hBAz1MsNIURhc+q05UjHXsB/tmo8LCyMFh07kjtnTlYuWpQ4gCQ1IIRg2qRJNGvfnq+//pxJk/7W24vB1zcjY8f+w48/9mTu3HN4eqbRW78+sdYVComMjubi7ducP3OG83fvcv7GDc4/eICzlHxbqRLlMmTQ6oWfAJ4DDfVapKI/UAlVdjVtU/cOKVSIk9HR5Jk7l3Sennzk58dHhQpRKlcuPipblkI5c2o9ITRDI5ok7AdD4ZHr18+yZMk01q8/887/kGRwJaDq+dO0Qyf69ehBp/b6hm44HkUKFeLQjh00btOGh48eMennn+2mx5U1kFJuBjYLIXyBX6WU/YUQrsAqKaVZWf+s2ktECJEH2JQQw9ZY/y1QHmgjdRgghOgH9APIlTNnuTsXL1rMrvvBwTRu04ZG9erx+48/ptqD5v6LeBo1qs6nn/Zi8OBh73yn7ZZ86NBueHunZ+jQP7XWlyBQ4eEvCAm5S2RkBN7eUURGRhIdHYWMfEZUdDSRUVFERUXxKjKSGxcvcv7qVe6GhFA4Z06VSObKxUdp0vCRnx9ZwzV6koSEvDuzeUQEvYDCwAhtBmmhMfCpeknoJZJovLqHSLyU3PHw4PyDB5yPiODc3bucDwnh/pMnFMmfn1KFC5OjQAG8vLzw8vTE09MT4ZkBLy8vPD298PLyIiLCk7RpffHzy4mbm7tO0X7z5jUDBlSkXbsv6Nmzx3vfawp2XFwcnTu3JmPGTCyaNsWiQ9btibCwMFp26kQOf38WzJhhcWfJIr1EChXS30skNpast25dBKI1Vs+SUs7SUt8kYKmU8rQQoi3gK6WcY45tKe5hCyG6o3oYWU+XWAOof/gsUHXrs1T7l69epXGbNgzp35+vhg61VLV2Sbp06Vi5cjP161chT568NGvWCtAdTh4yZCo9e5akcuWmZM2aiwcPgnjwIIiHD28TFhbE/ftBBAcH8ebNa/z8cuHrmxYvLy88PDzx8vLCx8NFJWoeHiqxi42ldf36fN++PYVy5sTVxeX98EO47v7HYcBaTEtAPABVXPtTPWWchCBvxozkzZiRlhqC/ioqiouvXnHuyhUexsTwMjych48eqS5CkZG8jHqTeHGKjIzkyZNnPHoUTMaMWciePQ85cuQlffq8ZMuWh2zZ8uLnl5fNm+eQOXMOGjfurtUWTS/7m2+GExn5isWLVxEh3CySjMoe8fX1ZfvatXTu3Ztm7duzeskSfHx0P2exY+4Z6tYnhPgF2CqlTHjS3BFVdMEsUlSwhRCNUTlLtaSUkSnZNsDR48dp1bkzE8aNo1vnzindfIqScAufM2culi1bT9u2TfD3z0HZsrodDx8fX776ahY//PAJvr5Z8fPLS758ecmePQ+VKpUnR4685MiRF1/fTO/NIKNVXHT1v9ZVLgmLgSaAKRPpNQMGA4GAqcM20nh6UtHTk4oJ6U61BKyThkbu3o3l0aNg9cXsNvfvB3Hlyl62bJnHw4e3kVIyY8Zxg97yjBl/sWfPDnbuPJzocYbjk2pF29PTk5WLFjHwyy+p26wZm1etIktm/cmyHA0hxEBUN32+QogCwCYgTEpp9p9qzW59y4DaQCYhxH1USa5GAe7ATvUBfFRKOcBaNmgSePEizTt0YOGMGXzcqFFKNGk3lC1bnr//nkPHji1YtOgAuXLp7nVQufLHbN6sypKnrSeEtoduesXaTC7Hx/M78K+J27mgiqN9AayVkmT16tYSZE4az86VywUXl9xkz577vc2N7fq3adM6/vjjF3bsOPRBTcDr4uLCzKlTGfPjj9Rt1ozje/fi5WV611J7RUo5HZieZHW/5NRpzV4inaSUflJKVyllDinlXCllASllTillafWSImIdHR1N5969mTB+/Ach1toekDVt2pIvvhjBF198QqzmTCwaaAq0sd3WrCHWO2NjqRUVxXighhnbjwIKAtWioriTkPfaXLT8FnN6jujan/fv32bQoL78999G8uR5v1ujvi6ZqQEhBONHj6Z0yZJ89e23tjbH7kmdT9uSMGrsWAoXLEiPLl1sbYpNadlyKOnTZ2D+/Ek6y+gTaoM9QjSy7pmERsUznj+na0wMqzw86GF6TQC4ouoK2NvFhSp373I0KsrwRvrcYTNFO2Ff6tqnUkrGjOlHz57DKVOmnGEbUzHTJk1iy86dbNq61dam2DWpXrB37tnDqvXrmTkl9T5110SfRyaEYNy4Wcyd+zu3bul+lGdMKESrWCeDuPh4ht28yZSwMA56elIzmYmNBPCFmxuzsmaleUgIyx8/TlZ92i5Gpoi2Ntasmc+LF6H06vVVcnefw5MuXToWz5pF36FDeZTc/yoVk6oF+9mzZ/QcOJD506eblTc6NZIjR14GDfqe777rTXy8cbMcWVusw1+/puXBgwS+esWRXLkoYMFuls28vdmdIwcjgoIYd+FC8gdsGBBtY3n0KIRJk0by009zDeavTu1hkQRqVK1Kr65d6TVw4AcxsMYcUu0UYVJK+n3+OR3atKF+nTq2NsfmaOpMly6D2Lr1P5YunUbXrkNMqseQWMfHxxP24gVPwsJ4GhZGzIMHSClVkyA8e6bKBSIl8U+fIoHYp08Zc/gwVdKl4+98+XC1gnf1kbs7x8qUoeW1a1zZsoVORYognj/HSQjE48cIIVTv799Xvb97lzSenmTOkIHMvr54J01mn+RhpK6cI7qQUvLDD5/RseMAihQpZfHf68iMHTWKqvXr88+cOQzs29fW5tgdqVaw5y9Zwo1bt/h37lzDhVMJxnpiTk5O/PjjXDp3rkadOs3JkSOPzrLacoQ8evyYeYsXc/fqVZ6EhvIkNJSnz5/zJDSUsJcv8UmThkzp05PJ1xcPVKEYJycnxOvXCHX7IiYGAYiYGPp/9BEDM2dW5cC2Etnc3AioW5eRN28y89o11UUEiHdze/ve1VV1MXF1JSIyMvF3xcbGksnXN1HAM2fIQOYMGShZsSKfduiAu7u7SaK9Zct/3Lt3kylTVryzXtvIxw8NV1dXls6ZQ7WGDalTsyZFC5s1gjv5GHrq/vo13LqVMrZokCoF+8bNm4wYM4aALVscIiWlJTD1tjlfvsL06vUVY8b0Ze7cHVrj+0nF+snTp0yYMoW5ixbxSaNGfFS4MJkzZCBT+vSJIpYhXbp3Z+HRNddW0j7aKRDE9XRxYWqdOu/+MH2TUaqJio5OvDAl3Dk8CQ1l9YoV/PDrr3wzfDi9unYFt3c31/aTQkOf8MsvXzB9+gbc3Iw/NhP+39TaL1uTQgUL8tOYMXTp04eju3enqrQRySXVCbaUkl6DBvHd//5H8aJFbW1OimBIrHVpYa9eX7F9+yoWL/6Tbt0+f+e7pJ7espUrGfzVV3Rs25bz69aRw9R0dQ6Mp4cHufz9yZVkpwzr0YOjZ88ydvZsfp86lfXLl5O7WCWd9cTFxfHdd31o3vxTPvqoos5y+kjNg2k06dujB5u3b+eXSZP4fpTNk3raDanuoePWHTt4FhrK4P79bW1KipCcB1IuLi5MmbKC2bN/Ze/eTVrLSCmZ/Ou3jPrhBwK2bGHaV199UGJtiMqlS7Nt7VrGf/cddZs14/CutYnfJb3oTZjwPyIiXjJs2M866zPmRuNDeAgphGDqb7/x54wZPH32zNbm2A2pSrCllIz+6SfGffut3cx3Z00sceLmzJmPv/9exzff9OTChVPvfBcTE0P//t3ZuHUrR3fvpmTx4slu7z1SImhr7QtMSAidP/mEtf/+S4/PPmPOnH/eK7Jkyd8cOLCVv/5ao9ziG0me3Lnp0KYNE6ZMsbUpdkOqComs27QJKSWtmze3tSkORalSlfjhh5kMHNiC5cuPUL58LkJDQ+nSpTVZMqRj39atKTtkOEkA+C6wj7czy+jCC1XiBm/rWaaXapUrc3D7dpq078DNm9f58cff8fd3ZunSTcyc+TPLlh0iXTpfi7T1oYRGvv3f//ioShW+HDyYbFmz2tocm+MQHnY8ToTj897yTpn4eMb89BPjvvkm1aZL1cTSt8UNG7ahR48v6d+/KWfOnKJevcqUK1eRVYsXvxVrUx8MJuNBYryU7ABaAaVRZc3ZZWCZC+QGhgKXzW7ZTNS/NX++fBzbtYPAwLN06dKGw4cP8M03Pfn777XkyKF9Rh0F3WT396dbp078Mund0bna9ECbLqQ2HNrD1vxzVq1ZjoeXDzUbtycc8UF4H5bm22+/JDT0FvXqVWbSpGn07NkPJ1370Uq9OkLj4ljw8iX/REaSBhgELMF4r/keqpy8dYGiUVEMDA+nZZYsuJp6ETf292kJ6fj6+rJmzTa++GIATZrUYvHiVZQqpfthpLmkdi874fwe/OUYKlQoxoCh35Ajh+kTH6cmrDqBgaUoW7a83LfvpM7vY2NjqVixOJMmTaNOnfp660oNB7ipXoQx2pOgO3Fxcdy+HUT+/AW0z81obIVJSZqrI0l9JwIDmX7+POuePKFZmjQMjI2lclQU5iYTeA2scXdnupMTN+Pi6OfnR9/SpfHPn9+4bn2moKO+l9KbW7dukj9/gcR1pvwXxpAajmcwfEyPHTuKsLBQpk7VOgncO6RNK5I/gUHNmvonMHj9mqxHj24zlA/b0ji0h53AkiXzyZbNj9q16xksq+3ASC0HvbloCoSzs/M7ApMcgp884f7Tp0RGRxP58CGRMTFEvn6ten3yhMg3b3j1+jW7L13iSXg4n2XNyrU8ecjs4pJsD94N6OjqSkd/fwJ9fPgnJITiW7ZQL08eiufJg5erK15ubnhlzoynm5vqfbZseHl44OXuTsEcOfBJZtxeCPHevjRmJKQpOKKXbU7YYujQ/1G2bCEGDRpGoUJFrGBVEgxdNaOj9X9vJRxesF+9esX48d+xZs02s5M7OZKIWzJGp++YTI53HRkdTc8JE9h1+jQF/P1VIggqUXR3V72+eYOXmxvpPDwYW6UKjd3ccH740Pgk0iZQMk0apteqxa9v3rDy5UvuC0FoZCT3X7wg8vnzxAtJFBAZE0NEVBS3Hz5kdNeuDE+YzMCC6BPt1DbS0VLHa4YMGfjqq2/57rv/sWLFRovU6Yg4vGCvXr2cChUqU6pUGYvW60gibghLe3X6CAsPp/m335I3WzaCV6zAI6ELm4GQSOLnbNmsItoAaV1d6V2ypFEhkbuPHvHxqFE8DA1lQv/++p0BzTHlGu91Td6b0Kyl/hN78bKt/cCvd+8BTJr0M0FBt8ibN59V27JXHL47xdy5M+jVK0XmQXDop9JJPTdreHIhT59S84svqFikCAtHjnwr1g5IrqxZ2T9lCocuXKDnhAnExsVZvI2U+E+sga3OA09PTzp16saCBe/Nc/vB4NCCffr0SUJDn1KvXkOb2pFSB7Cl6rSGMFy/f5/qn3/Op/XrM+mzz1JF18oMadOy8/ffeRwWRuvRo4m0UdzSVtijg9KzZz+WLJlPTEyMTe2wFQ59Vs2fP5MePfrZ7ahGeznIEzBWqE2NX5++do1aw4bxbZcujOjUKVVNFJHG05P1P/5Iem9vGo0YwfOICIvWn/CfJPciaqljzN6O2aQULFiYYsVKsHHjWsOFUyEOK9gvXrxg/fpVdO1q9ozxKYq9nASW9q73nD5N45Ejmf755/T++OPkVWbtmICZMQhXFxcWjhxJ+UKFqPnFFzwwlNtC48JmTGzZ1qEQexfppPTqNYB582bY2gyb4LCCvXz5YurWbUiWLI41XNXeTw5THl6t3r+fjj/+yIoxY2hVvboVrTIDb8sOUHdycuKPgQPpXLcu1YYO5fr9++8WcMA5vuz5ONRH06YtuX79Klevpvh4VpMQQuQTQswVQixXf/5TCDFbCLFSCGHWAeqQgi2lZN68lHvYaC1MiQ/a24k1a9Mmhvz1F9t/+43apUvb2pwUQQjByM6d+aZzZ2oNG8bpa9eM2s5eenDYUyw6Obi5udG1a2/mzTM8iMaWSClvSSl7a6zKKqXsC5wGzBqy6ZDd+o4fP8qbN2+oUaO2rU2xOLY8kYydq/HElSv8MH8++//8kwLZs5vXmJH92iRwBNjOu8mfhJb3XkBHIJd5FhlNn6ZN8fXxocU333Dz339xN6I3jL4ufpbCkUXYVHr06EuNGmUZP36CdbIfGopTRUYC5BRCaA7BniWl1NeF5YYQYi2qQ/Z3c8xySA97z54dNG3aMlU93LI1pniBq/bvp0+dOuaLtRG8kZJlQGWgGxCLyrtwAZxRHbhOvCvcQUAZVKJ9zArd8DRpW7MmOTNk4OCFC9oLOGCIxJHIlSs3uXPn4fTpE7Y0456UsrzGolOshRCZgExSytbASqCZOQ1azcMWQsxDZdRjKWUJ9boMwH9AHuA28ImUMszUug8d2sfgwcMtZ6yCSWw/cYLpn39uuKAmhgbEqD3u0DdvmBUayrTISAoA36A6iIztB/Qbqqx9HaOj8XN3Z1h4OK2zZn17oFvwCV+jqlXZfuIE9crqTTuRSEp42R8S1avX5sCBACpXrmZrU7QihEgP/AqUB3oDTkKI6UA2wLTZr9VY08NegCo9sSYjgd1SyoLAbvVnk4iJieHUqeNUqWJnD7kcGFO864ehodx59IiKFp5+7WpoKANPnqTAiRNcef2aDR4e7PX2piXGizVAWmAYcN3Li+G+vkwNC6PA8eNMunKFFxbuu9uoQgW2nTDNw7OHeHZqoUaN2hw8GGBrM3QipXwupRwgpSwgpfxNStlXSjlQStlGShlsTp1WE2wp5X4gNMnqlsBC9fuFqNIdm8SpU8cpVKgI6dKlS56BCoDpArLz5EnqlimDi4X6vp+8d4+ma9dSc8UKMrm7c6l8eRZky0aZZNbvIgRtfXw4mCsXK4sV41RoKHnnzOHzdet49uqVRWyvUKQI9588IeTpU5O2U0TbMlSpUoMTJ47y+vVrW5uSYqR0DDurlPIBgPo1i66CQoh+QoiTQoiTT58+SVx/8OA+qlWrZX1LPwD0CoeOGOz2kydpVKFCstuOj49nUkAAH8+ZQ8v8+bnTpw/jSpYkmxUeIFXw8eHfqlU5360b8VJSZvJkDrx8mex6XZydqV+2LDtO6k79q2A90qdPT/78BTlz5sPZ/3b70FFKOSshmJ8pU+bE9YcO7aN6dUWwbUF8fDw7T52iUXmzUw0D8PTlS5pPnMiqK1c4/vnn9PvoIzxcrN9hKYePD3+1bs3Mtm1pP3UqP65dS1x8fLLqbFShAttNDIuA4mVbioQ49odCSgv2IyGEH4D69bEpG79584YTJ45SpUoNqxj3IWGOYJy7eRNfb29yJ2NS232XL1Pmm28okTMn+8eMIU+GDNofBFrq4aAWW5sULcqpn35i14ULNPzlFx6EhZk9eUGjChXYeeoUcdp6pRjoKaKIdvJRxbH32dqMFCOlBXsD0F39vjuw3pSNz549TZ48+fD1tcxEph8q5grFnjNnjO4RoY3lly7R8a+/mNWnD7916oRrCnjVusieIQO7v/2WGoULU/bbb7n76JFZ9eTInJlM6dJx7uZNC1uoYAxVqtTg+PHD2i+YqRBrdutbBtQGMgkh7gPfo+riskII0RvVZNjtTanz4cMQcudWJjK1FZnTp+fYZfOHAy/dtYu/uneniZ2MjHR2cmJsu3bcevyYbSdO0K+Z6V1j4+PjefriBVnMdCKUrn7Jw9fXFxcXF168eEGGDBksV7GhO7xw29wdWbOXSCcppZ+U0lVKmUNKOVdK+UxKWU9KWVD9mrQXiV5evnyJj09aa5n8QWCUd63jVr7WRx8RcO4cRs8DqnHQSyk5evkyVQsVMm7bFKRa4cIcvnjRrG0v3bmDr48POTJnNlxYB0poJHl4e/sQEfFh7EO7feiojfBwRbCTQ3KFIXe2bHh7enL5zh2Tt70ZEoKnuzv+9hDOSuI9VSlYkCNJBdvIGHrA2bPUKlVKdwEjRzwqom0+Pj5pCQ9Pfq8fR8DhBDttWkWwTUEz3Y8lqF2qFAHnzhkumETwjly8SGULD7YxGgPiWzxHDh6GhfH0xQuTtgMIOHeO2voE2wQs+T99SCgetp2ieNjGY/bJb8AjrF26NAFnz+qvQ4vQHb18mSrFihlnQzJ6oZiDs78/FYsU4eilS6oVuuZ8TIKUkn3nzun3sM3A0hfZ1I7Kw/4w9pUi2KmIlDjRDcaxdUwScOTiRSoXK/a+GNs6e7+aqsWLmxzHvnT7Nmm9vMiZJYv+35GMRFCKcBvGx8dHCYnYI8pDR+1Y7KQ2QlhyZ8tGGg8Pk+LYr6KiuHr/PmULFkyOdValSrFiHEnwsJOiQ4wDrOBd60IRbt34+KRVQiL2SHj4S7wtPJOII5PiJ7FauGqXKsXOU6d0fp+UMzduUDx3bqPyRtuKysWKcfzKFeJ1jXzU8tt2njplvGBbKN2qItzvkyaNNy8tkGrAEXAowc6RIxe3b9+ytRmpExMEpW/TpkxcsYIYzaQ7ekIC6dKk4YWFEi4ZhRkX9RcREaT18kLoy/Gt8RsDb93i6KVLtKpmm9Seimi/JSjoJrly5bZspf7++pcUfs6SgEMJdqVKVTl27LCtzbALbHnCVi1RgoLZs7PxyBHVCgNx6GK5c/MgNJRnSXth2BFHL1+mUtGiRk+K8c+GDQxq1Yp03t7Gx+EtPKmBItqqgUvHjx+hYsUqtjYlRXA4wT569JDxAzcUiIuL49KVK8xfsoTPhg3ju3HjeP78+buFTBEStTj1bNyY+du2GbWJc86clC9UiONXrhhdf4qh9pSOXb5sXLdDf3+iYmL4LyCAHo0aJbv5a9ev03fIEEaMGcOaDRsIVmaqMYnr16+SPr0vWbPaxuNNaRxKsHPnzgPA3bumD9xITejzrIJDQli7cSOjxo6lbrNm+ObKRYsOHdi1dy9FChbk0ZMnFC5XjqnTpycrj3DbmjU5fPEiwUbmra5crFiyhrUnGwMXgqOXLlHJyH7i6y5epHyhQoZ7h2hDLcjPnj3j86+/pmqDBuTMnh3vNGmYu2gRpapWJUeRIrT99FMmTJnCvoMHeaUnnPShe9nHjh2mUqWqtjYjxXCoSXiFEIlhkQTx/tDQdoI+ePiQsb/8wubt24mOjqZS+fJULFeO/w0dSsVy5ciYMeM75T//7DO+Hj2av2bO5JfPP6ddo0amzY/p749XSAjtmzRh8fr1jOzXz+AmlYsWZfqGDaoPhqYLS2Fev3nDuVu3KF+4sFHl561ZQ5/Onc1qK+b1a/76809+mzyZDm3acPnkSTJnypT4vZSSW0FBHDt5kmMnTzJizBgCL12iQL58dOvUiSH9+1tn0lkH5UMTbIfysEGJY2sSHR3NzxMnUrJyZdKnS8e+LVt4EhTE5lWr+H7UKJo0bPieWAOUKFaMLatXM3PqVH6eOZNqnTpx+PRp0xr396dnmzbMX7PGqBBVpaJF9ffCsCHnbt4kv78/Pl5eBsveCQ7mzKVLtKxXzyTvWkrJf1u2UPTjjzmwZw8Htm/n70mT3hFrUDkl+fPlo/MnnzB1wgSO7tlD6J07/DN5Mnv376dEpUps3Lr1nX3+IXvZimDbORUrVuXEiSO2NsMmJJyYUkpWrVtH0QoVOHnmDMf27OG3cePIny+fSZ5yvUKFOLV6NQM6dqTj8OG0HTKE67dvG7195dKlEUJw+MwZg2WzZshAujRpuH7/vtH1A2b1+DCVY5cvU6lIEaMEeOG6dXT8+GM83N2Nrv/w6dNU7diRCXPmMPfHH1k/fTpFTEiC5e7uTtVKldi0ciV/TpjA16NH06hVKy5qhJg+RNF+9uwZDx+GUKxYCVubkmI4nGCXLl2Wa9euEBZm8mTrDk3CCRl0+zZ1mjZl/IQJzP37b9YsXUr+fPnMrtfJyYlurVpxdetWKpQsSZWOHfn+zz+N2lYIQa+2bZm/Zo3hwv7+VC5a1OyseCZjgvd7+OJFo+LX8fHxzF+zhl5t2xpVb8zr13QYNowOX37JZ506cWLVKupUrqz60syHi40bNOD8kSM0a9yYOk2b8sWIEcRYeHJhR+Hw4f2UK1cRZwvNL+oIOJxgu7u706lTN8aM+drWpqQ423ftonK9ejRv0oTTBw5Qt5blpkrz9PBgZL9+XNy0iflr13JQ28AYLXRt0YLVO3YQYUQ/61bVq7Ng+/bkmmpRwsLD2Xr8OE0ThFQPAcePk9bbmzJG5kSZsnAhL8LDubp1K91atcLJyTKnm6urK0M/+4zLJ05wLziYmo0bcz84+IPysqOjo/nhh2/o2bO/rU3RiRCilhDiPyHELCFEdUvUafAIEkJkEUK0FkIMEkL0EkJUFELYVOh/+OE39uzZwZ49O21pRooRHx/PT7//Ts+BA1m5aBHDhwyxmleRNVMmfv3yS774+Wej4s1+WbJQo1w5VmkT4iRebpsaNbj14AFnrl+3lLnJZs7mzTSrXBm/pLF+LR76vNWr6d22rVFhpwePH/P73Ln8PXo0Xp6eljL3HTJmzMiqxYtp06IFFevUYd/Bg1Zpxx759dcfKFq0OK1atbNOA9my6V+yZAHImTBRuHpJ+vS9HTAM+AwYagmzdAqvEKKOEGI7sBloAvgBxYDvgEAhxA9CCJsk9kibNi1Tp85i6NC+qT5L1+3bQXzasSmbtm3jREAANS01sk7PLXmnZs1wc3Vl4dq1RlXVq21b5q5e/e5KLYLn6uLC4FatmGpMCCUFiI2L4+916/iineGT/vnLl2wKCKBL8+ZG1f3tlCn0bteOArktPAIvCUIIRgwbxsIZM+jQowd//j6GqKgoq7Zpa06fPsnixfOYNGmarU25lzBRuHqZleT7P4ExwHiwzLRC+jzlj4G+UsoKUsp+UsrvpJRfSSlbAKWAM0ADSxhhDvXrN6JmzbqMGTPCViZYlWfPnjFq1JfUrlWOCmXLsm/rVrKn0KASIQRTv/2Wb6dMITwiwmD5prVqcf3OHa4FBRks27dpU9YfOsTDpIN3bMC6kyfJmSUL5XQ9ANTY38s2b6ZhtWpkNGIChpOBgWw7cIBvBwzQX9CCg2Qa1K3LsT17OHX2LOXLF+HffxelynkOY2JiGDiwJz///AdZsmS1tTl6kVJel1IOQDU1okmza+lCp2BLKf8npbyr47tYKeU6KeVqbd+nFD///Adbt25g//69tjTDokRHRzN16u+UL1+E+OiXXDpxgjEjR6Z439sKJUvSoGpVfpmV1Gl4H1dXV7q2aPH24aOeC0uGtGnpUKcOM3btspSpZjNl61a+MPIB4rzVq4162Cil5POff+bHL74gbQonKsudKxdrli5l2dw5zJ8/k5o1y7Fr1/ZUNTL4999/InfuvHzyiXn94FMSIUQFIcQMYAYqLzvZGBPDTi+EGCqE+EMI8WfCYonGk0v69OmZPHkGgwf30TsazBGIj49n2bLFlCtXmGPHDnNw+zamT55MVlWszCb8/OWXzFqxgiAjuuL1bNOGhevWERsba7Ds0NatmbF7NzFv3ljCTLM4eesW90JDaVXdwLMgf3/OX73Kw6dPaVDVcH/f5Zs3Ex0TQ4/WrS1kqelUq1yZIzu2MnLk93z99VBatWrIuXOGu17aO+fPn2XevBlMmTLDtIFeNkJKeUJKOUBK2VlKaUReBsMY8/BwC5AHCAROaSx2QZMmzahUqSrffjvcYT2JvXt3UbNmOebMmc6cOUvZ+O8ik/rpmoyRt+LZs2bli+7d+d+ECQbLFitQgNz+/mwzottesTx5KJ07N8uPGOhPb8UQ0NRt2xjcoAEu2h7eJml3/vbt9Gjd2uCD3sioKEZMmsSUb76xWI8QcxFC0Ll5fY4du0CLFm1p1+5j+vT5lDt3btvULnN59eoVn33Wg/HjJ+DnZx+TXtgCY44qDynll1LK+VLKhQmL1S0zgd9//4ujRw/x/fcjiY6OtrU5RhEfH8+ePTtp3boxw4Z9xldffcuuXYepUsUivX8sxvCePTl54QL7T5wwWLZ3u3bMXbzYqHq/aNKEyVu32mTkY3BoKJtOn6ZPnTpvV+q4OMTExLB0xQqjPOaJ8+ZRtUwZapQvbylTk00G12h69x7A6dPXyJevADVrlmPEiC+4efOGrU0zmocPH9C5cys++qgMnTt3t7U5NsUYwV4shOgrhPATQmRIWKxumQmkT5+e9et3cv36VSpWLMbatSvt1tt+/PgRkyf/Rpkyhfjuu69o3rw1x49fpFWrdggh7K4vraeHB6MHDmTCnDkGy3bo1YuAgwd5+OjR25U6hLBhyZJ4ubkxc/duS5lqNJ8vWsTABg3wNSLGvGHLFooXLUp+A/20o6Kj+WvJEsYPNbH3Vgpk5/MhHB8fH775ZizHj1/E1dWVBg2q0rJlA9auXZmsJGDWJCoqiokTf6Zy5ZKUKlWWv/6a7RChEGtijGC/Bn4HjvA2HHIyOY0KIYYJIS4KIS4IIZYJITySUx9A1qzZWLZsHX/9NYeJE3+iceOanDpl2CtMCeLj4wkI2E337p9Qrlxhrl+/ypw5Szl06Cy9er1N5mNvYp1A52bNOHb+PDfvan0GnYiPjw+tmzVj0bJlBut0cnJiTt++jFm1invW6jGiJcn8msBAAu/dY7SRMeY5CxfSp1s3g+WWbd5MxY8+omCePKZamSIkHFtZs2bjxx9/5/Lle3z6aS9mz55GsWK5GDt2FEFB9jE5iJSS1av/o0KFopw5c5Ldu48ybtxvuLg4VK46q2CMYH8JFJBS5pFS5lUvZo+FFkJkR9WJvLyUsgTgDHQ0t76k1KpVl/37T9GlSw86dWpJv37dCA42MX+FhXj69AlTp/5O2bKFGTnyC6pWrUlg4G2mT59HhQqVHMZb8PTwoGebNkz/91+DZft0787cRYuMusMpliMHI5o3p/Y//3DEyp6mlJK5gYH0X7WKuX374mFEr5s7d+9y8swZ2rRoYbDuv5YsYXCXLpYy1+q4u7vTvn0ntmwJYMuWAF6/fk3dupVo1aoR69ev5o2NHgifOnWCRo1qMHnyr/zzzwKWLl1D/vwFUt4QB55x5iIQaeF2XQBPIYQL4AVY9Gx1dnamW7fenDp1lRw5clK1ail++ul7q/ckefPmDceOHWHSpF9o1aoRZcoU5PLli8ycuYgjR87Tv/9g0qdPr3Vbe/WuExjYqRML163jVaSOQ0Ed+qhSUZXb4aChB4pqvmrWjEnNm9N6wwbGBgYSa4WY9tOYGNps2MBfZ88S8NlnVC9SxKjt5i9ZQuf27fFMGKmoI7xz+MwZXkVG0shQjxMbo+sYK1SoCD//PInLl+/RqVM3/vlnKkWL5qR//+4sWTI/RR5UhoQE069fNzp1aknXrr3Yt+8kNWrUtnq7joYx9xhxwFkhxF4gMcuMlNKsoZZSymAhxETgLhAF7JBS7khaTj3Msx9Azpy5zGkKHx8fxoz5iR49+vH99yMpXjw3FStWoUyZ8uqlXLJmqnjz5g1nz57mwIG9HDgQwPHjh8mbNz/Vq9emb9+BLFy4gnTp0hm2087FGiBPjhxUL1eOpRs30q9DB53lhBD06d6dOQsXUsOIbnAArUqUoJKbGz02bKD67t0syZcPS/lU20ND6XXiBJ2LFWN506a4G+kZxcXFMW/JEjb+95/Bsn8vWcKgLl1s3jMkuXh4eNChQxc6dOjCzZs3CAjYxa5d2xk7dhQeHh5Ur16bGjXqUKNG7WTPofjw4QNOnz7ByZPHOX36BKdPn6BPn884deoqPj4WGRSYKhGGbl2FEFofy5rbU0QI4QusBjoAz4GVwCop5RJd25QtW17u25essDkA9+/f4/TpE5w5c5LTp09y5sxJvLy8NAS8PB99VBonJydevnxJePjbRfPzy5cvOH/+LMePHyZXrjyJB3HVqjXJkMH057EpLthmhh92HT7MsF9+4fyGDe+HczS8zydPn1KwTBluBwaq7igS2tNsV3MCg5AQCAkhPjiYv65d48cLF/gtQwZ6pk2LePAAjBht+Q7e3kRly8bIyEjWPnvGgqpVqVumzFs7NT3lBAHXXOfvz7adOxn944+c2Lfv3bqT7LuQR48o0aIFQbt2kc5coUnhadHCTRwlLaXk2rUrHDgQwMGDARw4EECaNGmoXr02+fIVIF269KRNm4506dKTPv3b92nTpsPbWzWjueqcO8GpUyqBjoyMpGzZCpQtW4Hy5StSoUJlMmXKbJHflzatOCWlNLurjhDilLxypay+Mo+fPSNrtWrbpJRNzG3HHIzxsC9IKd/pdy2EMC6hgnbqA0FSyifqutYAVQGdgm0pcuTISY4cOWnRog2gOhBv3w7i7NlTnDlzkj///J3AwHM4OTnh45OWtGnT4uPz7pI2bVq8vX3o0aMvs2cv0TpBgCk4gnedQL0qVYiNi2P/iRPUqljx7RdJBCdzpkw0rFuXZatW8VmfPu+W03OxcBKCzwsXpq6TE10CA9n86hWzpMSkPeztzbm4OLrcvUvxdOk4V7YsvllNH8I8Z9Ei+nQ33IVs1ooVdPz4Y/PF2gb4EG6SaAshKFy4KIULF6VPn8+QUnL16mUOHtxHcPA9goPv8+LFc16+fMGLF88Tl5cvXxAdHY2HhwcffVSGsmUr0Lr1J4wf/zt585qWu11BhTGCPVsI0V1KGQgghOgEfAFsNLPNu0BlIYQXqpBIPZLZ68RchBDkzZuPvHnz0bp1+xRv35HEGlT7a3CXLvy1ZMm7gq2FPt26MeqHH94VbCMpmSYNx3Pl4punTyn14gXzMS5pTRww+fVrJrx5wx9Zs9KlYEGzROHxkyfs3rePedO0JBfSuOi8fv2aWStWsHPePJPbsDWmirYmQgiKFClGkSKG08y+efMGIYTSw8NCGLMX2wGrhBBdgOpAN6ChuQ1KKY8JIVYBp4FYVEmkDCesSGU4mlgn8GmLFoz64w8iXr3CO00aneXq16nDgGHDCDhwgNr58+uuUIfH7eHkxB9ZstAkOpre0dGUAjRTQ2oL5F0D0sTGcszTk7xp04KZHtzfM2fSulkz0qbVn4xy34kT5PL3p3jBgma1Y2uSI9rG4urqatX6PzQMCraU8pYQoiOwDrgHNJRSJit/o5Tye+D75NThyDiqWAOk8/GhcqlS7Dh0iDYNG+qMvzo5OTH1t9/o8dlnnFm5El8jHr5qo4GLC+eArbwv0knl2Ado6umJUzJutY+ePcvM+fM5mTR2rYX1u3fTsm5ds9tSUDAVnYIthAjk3XMkA6o+08eEEEgpP7K2cQr2SYu6dVm/e7dKsPXQvEkTdu3dS9/Ro1k5dapp4QmNmdV9vb3pbOyDx2SI9fOXL+k0fDgzp04lZ44cestKKdmwdy/bZs82uz17ICW8bAXLoc/DbpZiVnxAOLJ3nUDzOnUY+/ffxMXFYWjemwnjx1OlVi1m/vcfAzpabHyUxZFS0m/MGJrWqkWrZgYOfX9/zu3ejZurK0X1hXscBEW0tWCo546NQj36BPuZlFKvWyOE8DZURuEtqUGsAXJnz07WTJk4fekSFXLm1FvW3d2d5ZMmUa1zZ6qWKcNHdtqbYs7KlVwNCmLRokVGld9+8CAf16yZ/J4OKdylTxeKaDsG+nr6rxdCTBJC1BRCJD5dEkLkE0L0Vk8f1tj6JqYOUotYJ1CvcmV2GzmasVDevPwxciQdv/ySV3Y4fdXFoCC+mTyZ5X/8gYeHcWltdh85Qr0qVaxsWcqS2o7R1Ii+GWfqAbuB/sBFIcQLIcQzVP2lswHdpZSrUsZMxyY1ngj1qlRh99GjRpfv2rIl5UuU4HNtXeWMwdvbuMUYNHqmRMXE0GH8eH776iujwxsxMTEcOXuWWhUqmPNLFBTMRu9YWinlFillF3Xip3RSyoxSyqpSyp+klA/1bauQuqlVoQJHz583Kf/4tNGj2X/+PP9aO6WqCWGGYdOnUzJvXnq2aWP0NkeOH6dYgQKkN9DtzxFJjc5FasKxkx84AHZ5Alggbpo+bVqKFirE8VPGTz7k4+3NijFj+Pzvvzl+w/YJ9Kft2MG+c+eYMWyYSbHovfv3U6dSJStaZlvs8phVABTBtiqp/cCvXaMGe/fvN2mb0gUKMPerr2g9eTJ3njx598sUfAC39exZfly3js0//0w6EyfLDTh4kDpNUjSFhIICoAi2QjKoXb06AQcPmrxdi2rV+KppU5pNnMhLQyEVK4j4+ZAQus+YweovviCfifVHR0dz6uxZqhmYgcYo7KSHiDZSu7OREggh0ggh5gkhZqkHHyYbY2ZNf2+SPm3rFN7iQ7j9H/AWEIvqVapw4vRps+bR/KJJE2oULkyHxYuJjYtLti3G8uDlS5rPm8fUbt2oasZEx0dPnKBE0aJ4m+iVOyJ2fwzbP22AlVLKfkArS1RojIddXPODEMIZKGeJxlMjH8xB7u9P2rRpKV6kCMdOmp67SwjBn927I4Guy5YRHRtreRvhnQvT/fBwGs+eTZ9KlehkZK7upAQcOEDtGjXeqzu18sEcz0kIT3S7tC8ReAPkFEKc1Fj6JakmB6p0HgAWmZlDp2ALIUYJIcKBj4QQL9VLOPAYWG+JxlMbH8zBrSFURsWxdQibi7Mza7p3J15K6qxcySMrzgh0/MEDKi9bRpeyZfmufn2z6wk4eJDalphZxoHE3iHuGG3DPSlleY0laRK7+6hEGywUftbXD/sXKaUP8LuUMq168VF37RtlicZTEw55QFtANMyNYyfg5ebG8k8/pXGePFRatoxzYWHJtikpy+/coem6dUyvV4+v69Qxe3Si1vi1AwlvcnHIY9y2rAE+EUL8A2ywRIXGZOsbpZ44N7dmeSmlad0DUjEf1IGcRKCqVa7MyTNniIqKejv3oYkIIfi+ShWKZMhAg127mFegAM2SOTEEqPKD/BAYyIKgIHa3a8dHmZM3o8mxkycpVrjwBz2FlTKE3XiklK+AXpas06BgCyF+RTWr+SVUOeJBlcVPEWw+MLHWQlp1f+xTZ89SPZlDtTsULkze169pHhDAiqJFqZVM2367d4+1z59zrGFDsiZTrAEOHD6sfZ5KAzPpaC3vwCiibTuMiau0BgpLKT+WUjZXLy2sbZgjkCrE2ljxSDoXogaVK1TgyPHjFjGnYsaMLKtenU8uX+ZiTIzhDXSwWAhmPHjA1lq1yGpkfhBDHD1xgiq6ZtoxZT8qKJiJMYJ9C1CmjUhCqhDrBAyJiIHvq1SsyNETJ95bb64XVjdrVibnz8/HwcEEv3lj8va7Xr3iq1u32FKiBP5mhmmSIqXUL9ig96KW+H0qIVUd/w6EvgkM/kIV+ogEzgohdgOJLo+Ucqj1zVNIMXTd1hshMpUrVOCr775DSmmxiVU7Z8nCvZgYPg4O5oCrK2mNrPdshgx0Dg5mVfHiFNMzhZkhwvF5R5Su37hBmjRp8PfzM7yxtn2ZisQ6ASU0kvLoi2EndK49hYWecKYWUq13oSk0JghM3jx5iIuL4979++QykB/bFL7OkYO7oaG0iYhgi4cHbgZE+058PM1CQvg7SxZqpk9vMTtAFQ6pbEp2PjP3paOhiHbKolOwpZQLU9IQRyDVCrUmZoiLECIxjm1JwRZC8GeWLLSNjaVaVBQ5DAj2mfh4hmfKxCdW6MVx5MQJqpiaTjUVC7UmqVG0DT1DfvYsZexIijG9RJLO7QjwApUH/qOU0kampywfhFgng4Q4doe2bS1ar7MQLPfzY8f9+xgawD4EqOvra9H2Ezhy/DjdO3WySt2pgdQo2vaIQcFGNWF1HPCv+nNHVBNWvwAWAM2tYpkdoYi1YSpXqMDI779P/GzJk9fDyYkWuXJZrD5TCQ8P5/rNm5QpVcpmNjgCimhbH2MEu5qUsprG50AhxCEpZTUhxKfWMkzBsShUoAC37961tRkWJeHB4/2QELL7+eHu7m5rk+weRbStizHd+ryFEInZ2oUQFYGEVGVmZewRQqQXQqwSQlwRQlwWQtjt5HiKd20ccXFxuLgYc/13PFycnYmPt0jung8C5ZyxHsacYX2AeUIIb1ShkJdAH/XEvL+Y2e5UYJuUsp0Qwg3wMrMeBTshNjbWfgQ7WzbVa8JDv2Q+/HNxcSFOEWwFO8CYXCIngJJCiHSAkFI+1/h6hakNCiHSAjWBHur6XwOvTa0nJVA8BeOJjY3FxdnZ1mbox0zhdnFxIdZa6V9TKUpoxDroGzjzqZRyiRDiyyTrAZBS/mFmm/mAJ8B8IUQpVP28P1cnStFspx/QDyBnTts9cFIwjrj4eMt62An9mLNlg4cWnu/ZROFWBFvBXtAXw04YJuajYzEXF6As8I+UsgzwChiZtJCUclZCntlMmZKfuMdUFO/aNDRDIsnyrOyw77Ii2OahnEOWR9/AmZnq1x8s3OZ94L6U8pj68yq0CLYtUQ4007GrGLaxJMS69RCODy7OYSk6jZmC7TF0U2eFtO1GYcycjoWEELuFEBfUnz8SQnxnboNSyofAPSFEYfWqeqhStyo4MMkWbFt61gbadnFxIU4RbLNQnB/LYky3vtnAKOANgJTyPKrBM8lhCLBUCHEeKA38nMz6FGzM6g0bKFyggK3NMMpr1okO4fby8iJ9unRs37XL/LoVFCyAMS6Rl5TyeJIsbMkK6EkpzwLlk1OHtVA8AtM5cPgwC5Yu5eyhQ7Y25S3J8NiTxuCdnZ2ZP3063QcM4NyhQ2S0wGw4CgrmYIyH/VQIkR91PhEhRDvggVWtUnAYXrx4Qbf+/Zn1559kzZLF9AqMEdbkeM0Wol7t2nzSujUDhg1DyqSpdRQUUgZjPOxBwCygiBAiGAgClCHpCgAM/fprGtatS/MmTUzb0JjuepopSs0VbQvGxn/+/nvK16rFkuXL6aokgjIapU+25TBm4MwtoL56ZKOTlFKJGSgAsGrdOg4fO8YZjVnTDZ6Yps5/aEd4eHiwdM4cGrRsSc1q1chtw4RUCh8mxqRXdQfaAnkAF42BM+OsapmCXRPy4AGDhg9nw/LleHt7G97AFKG2Y1EvVbIkXw0dSvcBA9i9cSPO9j66U8GuEEJ0QJXhNBqYKKW8Ysr2xsSw1wMtUT1ofKWxpDqUB47GIaWk96BBDOzTh0qmJvW3FRZ8CDl8yBDi4+OZMn16cq1S+PBojSotx9fAcFM3NiaGnUNK2djUihVSL/sOHuRmUBAbVxhOJZN0bsQUwcp9up2dnZn7999UbdCAfj164GOFGW5SG44WxzZ0g/fiBQA5hRAnNVbPklLOSvgghCjJ+wnyxgF/Aw+BDKbaZYxgHxZClJRSBppauULq5JdJk/j6iy8sP7JR1+S1dhgeKVigAPVr12bm/Pl8NVSZj/oD5Z6UUufTdrVmNtPy1XEhRAFgsKkN6kv+lDA1mAvQUwhxC9Ws6UJli/zI1MYUHJ9TZ85w8coVNnRM7tipFMDKnvbIL7/k43btGNyvHx4eHlZtSyF1IIT4GFUM2xv4n6nb63ORtF0ZFD5wfpsyhS8HD7b+7CvW8qwtKOKlSpakdMmSLPz3X/r36mWxehVSL1LKLcAWc7fX+dBRSnlH32JugwqOy7Xr19m7fz/9evSwtSm6SeGcJKOGD2fClClKNj+FFMGYXiIKCgBMmDqVQX37GteNzxQMiWxyRdiKIl69ShX8/fxYuXat1dpQUEjAwfJhKtiK+8HBrNmwgetnzli3IV2hEFNF11B5Cw53H/Xll4wcO5aO7dqRJOeOgoJFUTxsBaP44++/6d65c8omPrLDyQy00aRhQ5ycnNi8bZutTVFI5SiCrWCQZ8+esWDpUoYPGWJrUyyPBS4KQghGffklP0+apCSGUrAqimArGGT56tV83LAhObJnt2zFdpCFz1K0a9WK+yEhXLpi0khjBTslJET/YulpRo1FEWwFg+wKCKBJgwa2NsN4NL1mXe/1bWMGzs7ONKxblz379iWrHgUFfSiCraCXuLg49h08SN1atSxXqTZxNEdk7Yy6NWuyZ/9+W5uhkIpRBFtBL2fOncMvWzb8UlH4wlrUrVWLfQcPKvM/KlgNRbAV9LJn/37q1qxpWyNM8bJt6J1ny5oVv2zZOHv+fIq2q/DhoPTDVtDL7oAAPuvTxzKVJScndgoNnkluRrm6NWuye98+ypUpk6x6FBS0oXjYCjqJiYnh8PHj1K5ePWUadKB4tS7q1qqlPHhUsBqKYCvo5NjJkxQpWJD06dPb2hTT0Sf+VozH165encPHj/P69WurtaHw4aIItoJOTpw+TdVKlWxngAN63L6+vuTMnp3LV6/a2hSFVIjNYthCCGfgJBAspVRSudohERERpE+XzrqNGDN7enJJYeFPny4dEa9S5Sx6HwyGHrVERqaMHUmxpYf9OXDZhu0rGCAyKgovL69k12PSg7yk4mqO2Bq7jZWE3MvLi0hbndEKqRqbCLYQIgfQFJhji/YVjCMyMhJPR59JxQZhlTReXrxSBFvBCtjKw56CatbgeF0FhBD9hBAnhRAnnz59kmKGWYrXr18TExNjazOShaU8bK2YMkzc3H7Y2kiBAUBenp4O72HHx+s8NRVsSIrHsIUQzYDHUspTQojausqpZx+eBVC2bHmbp0CTUvL4yRNuBQVx6/ZtHj5+TGhYGM9CQwkNC3vvfXR0NM7Ozvhly0ah/PkpVKAABdWvhQoUIHeuXDg7O9v6Z+klMjISL09P6zdkTBzbgR5Aenl5ERkVZWszDJJwTF+5do2r16+/83rn3j1cXV3J4OtLxgwZyODrm7hofs7h70+BfPnInSsXbm5utv5JqR5bPHSsBrRQT0bpAaQVQiyRUn5qA1veIT4+npu3bnH95k1u3b793uLh7k6+PHnImycP2f38yODry0fFiycewJqv3t7exMXFcfvOHa7duMG1Gze4ev06G7dt49qNGzx6/Ji8uXNToWxZfhs3zi6Hfuv1sENCrCuilpjT0doir2Mf2LOHff3GDX6dPJmLly9z5fp1XJydKVKoEIULFqRIoULUrlGDIgULkjdPHt68eaPVGXkWGsqz0FCu3bjB3fv3uRkURHBICP5+fuTPm/edpUC+fBTMn996d2ofGCku2FLKUcAoALWH/ZWtxPrBgxBOnTqeuJw5c5L06dJRpGBB8uXJQ768ealepYpKpHPnJp2JPSZcXFwokD8/BfLn5+NGjd75LioqiptBQfy3ejVlqlfnn8mTad28uSV/XrKJjIpKGQ87pdB1UbSwsNujhy2lZNb8+Xw3fjxfDh5Mr65dKVywIJn0TEjh6uqKl5eXUWl137x5w527d7kZFJS4HD52jBu3bnHr9m0KFShAmXKVKVu2AuXKVaRo0eK4uCgDrU3lg9pjz58/5+DBAPbs2UlAwC5CQ59RvnwlypatwKBBX1KrXFEyZ8qUIrZ4enpSolgxShQrxscNG9K5Tx8eP3lid7Nvp1giI21hkeR42TYMoUgp7S4GPGrsWLbt2sX+bdsoWriwxet3dXVNdE6SEhMTw7nAQA6cvsSRIwf5669JPHnyiOrVa1OnTgNq165PgQIFP4jp1YQQFYHhQJCUcqQQwglYBEQBrkAvKaXOg8emgi2lDAACrNlGXFwcu3ZtZ+HC2ezbt5sKFapQt24DFiz4jxIlPsLJSfXc1Ydwa5qhlyqVKrFz3TqqN2pE/rx5qV+njs1s0aRE0aKcv3iRJg0b2toU00gJsdZzITkXGMjAvn2tb4ORzFm4kNXr13Nk9269HrW1cHd3p2L58hQtXwcYBMCjRw8JCNhNQMAuJk36GTc3dz79tCefftoTf38LT5RhBob8hOhoAHIKIU5qrJ6lfvamEynlcSHECGCAepUXECOl7CuEmAakAd1ilGo97Dt3brN48TyWLp2Pn192unfvw8yZi/DxSV5yH2tRIH9+/luwgE+6d2ff1q0UKVTI1iZRrnRpNlp6nkJ9XrOlveykdSet1wpIKTl97hxlS5WySv2msjsggG/HjePAtm02EWtdZM2ajQ4dutChQxeklJw5c4pFi+ZQuXJJqlSpTo8efWnQoIm9h03uSSmb6PpSCFES+CXJ6qTh30hVUbEaCJVS6vUcU9XQ9NevX7N27UpatWpErVrlefnyBatWbWHPnqN0797HbsU6gVrVq/PrDz/QrH17nj57ZmtzKFemDKe0zZJuCQEF7aKpLc6c3G59yXyga8rd1/3gYIQQZFfbkdzsf8nh8tWrdOrVixULF1KoYEGb2ZGArv0ohKBs2fJMmTKDS5fu0qxZKyZO/JkSJfIwfvxo7ty5nbKGWggpZaCUslmS5XmSYmWB61LKtsBtIURpfXXa9eXLWIKCbjFnznSWL19M0aLF6d69L8uXr8fDAQd99Pz0U65ev06bLl3YuX497u7uKdp+OD6JJ1ahAgV49OQJz58/154Ayto9RTQxt50U7n1z6uxZypUubfN47OMnT2javj2/jx9PrZTKtmgBvL296dq1F1279uLixUAWLpxDrVrlKV26HD179qNp05b27nXrRQiRHxgHFBNCXAeWAV8KIaYDmVCNUdGJQ3vYFy8G0rt3F+rUqYizszM7dhxi06Y9tG/fySHFOoGfv/+eTBkz0v/zz206C7ezszOlSpTg9LlzRm9jkWcBdtjFUSdJ7jZOnTlDudKl31mX0l52dHQ0rTp1onO7dnTv0iVF27YkxYuXZMKEqVy5cp9OnboxbdpkypYtxKxZ0+y226QhpJQ3pZRdpJRlpJRzpZSRUsrOUsqBUspPpJR6k9A4pGAfO3aYTz5pTqtWDSlR4iPOn7/FuHG/kT9/AbPqs+UDR204OTmxeNYsAi9d4rfJk1O8fU2BKVe6NKfOnk1xGywi2jYQ/lNnz9p08oL4+Hh6DBhArhw5GPfddzazw5J4eHjQoUMXduw4yOzZS9m7dyclS+bl11/H8cwOQocpicMItpSSnTu30aRJLfr2/ZSGDT8mMDCIYcNGkDZtWlubZ3HSpEnDhuXLmTZ7NqvXr0+RNrV5guVKl9Yex7Z3bCDWUkpOnT2r9YFjSnnZ3//0E3fv32fBjBmJPaDsieQ6R5UqVWHZsnVs3bqPe/fuUKZMQb7++nPu3r1jIQvtG/v7R7Xw/HkYtWqVZ/To/9GjRz9On75Gnz6fOXTYwxiy+/uzftky+g0dSlhYmE1sqFqpErsCAjh09KjZdZh9klpLdK0Ud58+ezYeHh7kzJHDKvUb4sDhwyz+7z/WL1+e6s+NQoWKMG3aXI4du4C7uzs1apSlX79uXL580damWRWHEOwnTx4zcuT3HD58jg4dulj0oYO9hUOSUrZ0aSqVL89uG007VSB/fhbPmkXrzp1Zsny55RuwRrImCwq9McdHbGwsQ776immzZ7Nn40abPXBcvX49/Xr0SLHBX/aAn58/48dP4Ny5mxQqVITvvvufrU2yKg7xuLVgwcJ8/HELW5thMxrUqcPOvXtp16qVTdpv0rAhezZtonmHDly5do1xvXrpvd32IdzoEMCuw4eJiIykZdGiuoXOQgL8KiqKhceOUbdyZYrky6e1jCm2A7wID6fDoEFIKTmya5fJ6Qssyc69e1nwzz82a99YTN3HxpA+fXq++uobANKmTf4F01AuMlvNAOcQHra1sHfvOoH6asG2JSWKFePYnj0EHDxIh2HDkp0rIzwign6jR9Pnu+/4Ydo0qg0fzpGL1rmdjYuLY96RIxTq2ZONe/dSvXNnpi5alOzh47fu3aNKx44UzJ+fzatW2VSs7wcH8/DRI8om6aGikLr4oAXbUShRrBhRUVHcCgqyqR1ZMmdm98aNeHp4UL1LFwKOHTOr2+H+Eyco1aoVcfHxnN+wgVOrVzOgY0c++eknPpkwgZvBwRazeceJE5QZNIj5a9aw5q+/2Dp7NkeWL2fF1q3U79mTO7GxOrfVdUF//fo1C9eupVrnzgzu0oW/Jk60ed/gXQEB1K1Vy+5T9ibgKM6SvaEItgMghKB+7dop5mXrO5nc3d1Z+OuvDOnShf7ff0+VevVYu3Hje95qQh2adUVHRzP8m2/o+OWXTP3mG+b+9BNpvb1xcnKiW6tWXN26lVJFilBpyBC+nD6d0Jcvzf4Ngbdu0XjECAZNn84PQ4awf8kSKql7bxTMk4f9S5bQqHp1yteqxbzFi4268ERERDBl2jQKlCzJ4g0bWDV1KgM7dzbbRkuyc88eGtata2szFKzMByvYjnaFb2AHYZEEhBD0bNuWS5s387/PP+fniRMpXrEi85cs4bWO4N7J06cpW6MG94KDOb9hA821iIuXpyffDhjAxU2biHR2pkiPHkxeuZIYEwKGD549o+/EidQfMYKmDRtyceNGWjdo8F583NnZmRF9+7Jn0yb+nDGDFh068PDRI611Pn32jO9/+om8JUty6Ngx1vz5J7vmz6da2bJG22VN4uPj2RUQQAM7SRpmLI52DtoDDvHQUQHq1a7NFyNHEhcXZze3vc7OzrRt2ZI2LVqwd/9+fps8mTE//cSwQYPo2707+Pjw5s0bfvr9d/6ZO5cpv/5Kx3btVOKpJx9J1kyZmPHDDwzt2pX//f47f69bx5DWrfHy8EBKmegNJ76qt7v76BFztmyh9yefcHXrVtIb0T+/ZPHiHN+7l3G//krpatX4e+LExIe7d+7eZdJff7Hkv/9o16oVh3bsUOXksFQuFQsRePEiaX18yJM7t61NUbAyH6RgO+KVPbu/P/5+fpw4dYrKFSva2px3EEJQt1Yt6taqxemzZ5kwZQo/TphAnZo1ORcYSMH8+Tlz8CD+fn4m1VusQAE2z5zJ7iNHWLF1K1LKRE9ZCPHe+zRp0nByzRrymNgP2s3NjR/HjKF5kyZ069+fKdOn4+bmxrnAQPp0787F48fNnhEooUeENY+57bt3O5x3rWAeH6RgOyotmzZl3ebNVhNsS4hK2dKlWb5gAQ8fPWL7rl2M+OILKpQrl6y+yfWqVKFelSrJts0QlSpU4MKxY2zduRMpJQ3q1LHI1FbWdhDWbtzID998Y9U2rIU1uvilZj7YGLYj0rpZM1avX2/ThFDGki1rVrp36ULF8uVtnrlOJ1oG7bi6utLi449p2bSpQ8xDGBwSwtXr16ldo4atTVFIARQP24EoV6YMMa9fc/HyZUoUK2Zrc5KHpSYm+MBZt2kTTRs1UmYstzCGDk09vUGtygfnYTti/DoBIQRtmjdn7caNtjZFwU5Yu2kTbVp8uKOAPzQ+OMF2dNq0aMGa1CLYNpwo16ZtW4hnz55x4vRpGtWrZ2tTkoUjO1EpjSLYDka1ypUJDgmx+ajHDxo7EfuN27ZRr1Yth4i1K1gGJYbtYDg7O9OyaVPWbtrE8CFDbG2O9bCUKGoLRtqJ4CaXNRs28Enr1rY2QyEF+aA87NRy69W2RQsWL1+e7ORFdoG/v/bFmvWnAh48fMj+w4dp1rixrU1RSEFSXLCFEDmFEHuFEJeFEBeFEJ+ntA2OTsN69XB3c7NOfuoUIjIykvWbN9Nr4EB6fvYZG7ZsISqZGQD18ezZMxYsXcon3bvzzQ8/cPzkSctc8Gx0ARg1diz9e/bUPjmyQqrFFh52LDBcSlkUqAwMEkI4eB+1lMXJyYmpv/3GN+PGERERkfIG+Ptz484dQk2c1eTJ06fMX7KEVp06ka1gQf6cMYPSJUtStlQpJk+bhl+hQnTo0YP/Vq8mPDz5d0PBISFMmzWLes2bk69UKTZs2UKjevWQUtJ9wAByFCnCZ8OGsW3nTmJiYoyuV0rJo6dPeWyj+QSPnTjBzr17+e5/qTtZf2pECNFdCDFbCLFBCFFSCJFRCDFDvdwVQujNpyBsPQhDCLEe+FtKuVNXmbJly8t9+04mu63UEhJJoEvv3uTPmzfFJluNi4tj3aZNTPzzT4KCgoiMicHF2ZmC+fNTIF++914zZszIzVu3WL9lC+s2beLchQs0qFOHVk2b8nHDhmTIkOGd+p88fcr6zZtZs2EDB48epVa1arRt2ZLmjRuTMWNGo2y8cfMmazZuZO3GjVy9fp2mjRrRpkULGtWr997DuavXr7N+82bWb97MxStXaFyvHi3VtqVNm5ZHjx9z/eZNbty6xY1bt9557+bqSlxsLB83asTwIUNSLA91fHw8VerVY2CfPg49I3pSLDnaMW1acUpKWd7c7YUQpwoVknoze8XGPubWrazbpJRNzGyjDNBYSvmL+nMmYJKUsrve7Wwp2EKIPMB+oISUUmcuTUWwtXPv/n1KV6vG6QMHyJ0rl9XaefXqFQuWLmXytGlkzpSJ/33+OS2bNsXJyYknT5+qxOzGjbeiFhTE9Zs3kVLi5elJ8yZNaNW0KfVq1zZ6rsHnz5+zeft21mzcyK6AAPLnzYuLgaRXL16+5MXLl7Rq1ow2zZtTu0YNoweUPHr8mI1bt7Ju0yb2Hz6MlBIPd3etF6IC+fLh6+vL8+fPmb1wIX/OmEGBfPkYPngwHzdqZNXJbxcuXcr0OXM4snu3XU6yay72Jtje3voFOz7+MZGRWS8C0RqrZ0kpZ2nUUxL4JcmmnwIRwGxgjJTynrrscOCclHKXXttsJdhCCG9gH/CTlHKNlu/7Af0AcubMVe7ixeTPipzaBBtg7M8/c/X6dZbNn2/xuh89fszfM2cyc/58qlepwldDh1K1UiWjtpVSEhoaSvr06ZOdXTAyMpILly4ZLOfh4UHxokWT3d6rV6948+aN0fHhN2/esGLNGib9/TdRUVEMGzSIrh074unpmSw7khIeHk7hcuVYu3QplSpUsGjdtsZBBdtkD1sI4Qr8BfwjpTynsX4z0EwaEGSbCLba6E3AdinlH4bKKx62bl69ekWR8uX5b8ECo8XUEJevXuWPv/9m1fr1dGrblmGDBlGwQAGL1J2akVKyd/9+Jv31FyfPnGFgnz4M7NvXYpPijhgzhkePH7NgxgyL1GdPfECCPQHVs7tLwG4p5UohRFWgnpRyvKHtU7wftlBlApoLXDZGrBX0kyZNGn75/nu+GDGCo3v2mH2bHBERwer161nw779cunKFQX37cu30abuZgdvYE9qWF2XNNLOXrlxh8rRpFCxThsb16tGra1fq1a5ttvd//cYN5i5aRODRoxa2WiElkVJ+rWXdYeCwMdvbIghWDegK1BVCnFUvH9vAjlRD508+QQjBspUrzdr+6vXr5CxWjNUbNjCkf3/uXrrEmJEj7UKsw/Exyfsytby1KFakCLP/+otb585Ro2pVvhk3jlJVq5pd39Cvv2bEsGFm5+VWSB2kuIctpTwI2CTfprUTydsKJycnxowYwdhffqFLhw4mb3/z1i2qVKjAhv/+s4J15pMc4bWX/zpDhgwM6tePQf364Zw+PbGxsSZP2Hvl2jVOnzvHegfud69gGVLPY+YPnMYNGvAsLIwTp06ZvG10TIzRvTdSCkt4yfbgaWvi5uamc85LfcyaP59eXbsqKVQVFMFOLTg7OzOgVy+mzZ5t8rbR0dF4uLtbwSrzsKTQ2pNou7u7mzRAB1T/zeLly1VzZCp88CjJn1IRvbp2pWCZMjx79szogSYAMa9f424ngm1PAmtp3N3cTJoBHmDVunWULVWKfHnzWskqBS0YEbK1TXdoxcNORWTKmJEWTZowb8kSk7aLjo62u5CIJbGXi4Cbm5vJHvbM+fPp36uXlSxS0E77MhER+4mIQOcSGbkOGJ3imbcUwU5lDOrbl3/mziUuLs7obewlJGJNYbUH0XY3MYZ98fJlbgYF0byJWaOfFcxmZSEYBehKDhYOTAfGp0s5m1Qogp3KqFCuHBl8fdm+S+8I13eIjomxeUgkJQTV1qLt7u5uUkhk1vz59Pr0U1xdXa1olUJSpJTXoSywQkeJiUB/9KXTsBZKDDuVIYSgdbNmLFu1io8bNTJqm/j4eDZt24Z3mjRUrVSJiuXK4e3tbWVL35KSQprS3f1CHjzg8LFjHDl+nPshISaldF25bh3r/v3XitYp6ObvLHDgMbQGNJ2ZEFSDtE+7wcAUt0rxsFMZl65cYeo//9C3Rw+jtxk6YADjv/uOFy9fMvrHH/EvXJim7doxd9Einlo5hagtvF5rt3n56lV+nDCBsjVq8FGVKixatoxMGTOy8b//KFm8uNH1DOrbl0HDhxMZGWlFaxW0IaV8Ah1QhT40GQuMRkr5JuWtsoP0qsZgqVwikDrziSTw5OlTKtety5gRI5KVevPly5eJmfJ27NlDudKladO8Oa2bNye7hRP22ypMYcnjQErJmXPnWL1hA2s2bCDi1SvaNG9OmxYtqF6litnD0aWUdO/fn8ioKFYsXJiqsvMlxZ5yiSQghPCCj15BAOALXAQGAfucDCVpshaKYKcSYmJiqN+iBTWqVuXn77+3WL1RUVFs372bNRs2sGnbNgoXLEibFi1o1bQp+fPlS5aI2DqmnJxjISYmhpNnzrBmwwbWbNyIs5MTbVu2pG2LFpQvW9Zi4mqt/9XesEfBBhBigVQJ9e9Ac2BTFSmlzRK6KIKdShg2ciT3goOt6om9fv2agAMHWLNxI5u2bSM0LIx8efJQIF8+8ufNm5grOn/evOTKmdPgEGxDJ2l8fDyBgefYt283Fy8GYuhYdXNzo0KFytSqVY88eQz3WzZ0LERERHAzKIibQUHcuHUr8fXGrVs8fPSIIoUK0bpZM9q2bEmJYsVQ5TWzPE+fPaNSnTpM+e23VNtjxH4FWzhDxVgYAGxFyhU2SauRaI8i2I7Pmzdv8CtYkDMHD5IzR44UazdB0JKK2c2gIB49fkyuHDkoWrgwJYoVo0TRopQoVozCBQvi5uam9QSVUnLjxnX27dvNvn27OXAggEyZMlOzZl3KlClHeLj+C4Cb2yuOHDnIvn278fLyombNutSqVY+aNeuQNav2pEk+hBMWFsbFK1cIvHiRC5cuceHyZa5ev87L8HDy5cnz3sWoQL58Rl2QLMmCpUvZuHUrq03sY+8o2KtgAwixQ6q869cFpJQ3LVWvWbYogu347A4IYOT333Ni3z5bm5JIdHQ0QXfucOnKlUQRvHDpErfv3iVfnjwUKfYRRYuWoFixEoSHh7Nv3272798DQK1a9dRLXfz9swMQEmK4zYTwupSSK1cusW/fHvbt282hQ/vw88tOzZp1qV69FhER4Vy6dIFLly5w5XIgL16+TLygaF5Y/P387CZu/OzZM/KVKsWDa9fem+osNWDPgg0ghMgmpXxoyTrNskMRbMdn8PDhZPf3Z9Tw4bY2xSDR0dGcvh7CpUsXuHxZJZru7u6JIl2gQEGtoQVjBBu0T2IeGxvLuXNn2LdvN4cPHyB9el+KFStB0aIlqFgsD7ly5rQbYdZH/RYtGNS3L62bN7e1KRbH3gXbXlAE28GRUpKzaFF2bdhAkUKFbG2OUZhzciZHsA3hKMfE9NmzOXriBItmzTJc2MFQBNs47N+tsCCOcmKawpVr13BxcaFwwYK2NsUorNkzxMI9Du2OxvXrsysgwODDV4XUizLS0cHZs28fdWrUsFoPBXshtYuxMeTNkwdXV1euXr/uMHdTCpblg/KwUyN7Dxygbs2atjZDIQUQQlC3Zk322NHDZYWURfGwHZj4+HgCDhxgyq+/2toUq/H69Wu2bNnA/v17DYYCPD09admyLRUrVkm1dxx1a9Zkw9atDOzb19amKNgARbAdmMCLF8ng60uO7NltbYrFuXbtCgsXzmH58sUUKVKMjz9uaXCKrNDQZ3z2WU9cXFzo1q0PnTp1M2oiB3uZ/9EY6tSsybBRo4iPj3eIni0KlkURbAdm7/791ElF4ZDIyEjWrl3JokVzuHXrBl269GD79oMUKGD8A9Wvv/6OQ4f2s3DhHH777Qfq129Cjx59qVGjdqoQuBzZs5PB15fAixcpVbKkrc1RSGEUwXZg9uzfT5dPPrG1Gcnm7NnTLFw4hzVr/qNixSoMGTKcRo2ampUHWghB9eq1qF69FmFhYfz33xJGjvyCyMhIunXrTZcuPciWzc8KvyLlqFurFnv371cE+wPE8V2OD5hXkZHUrlHD1maYzYEDAZQokYdPP21Dtmx+HDp0lpUrN9GsWSuLJO339fVlwIAhHD58jrlz/+X27VtUrFiMWrXKExsbm+z6bUWjevWsnvY2pXGUkJStcYiBM0KIJ8AdW9thBJmAp7Y2wgwUu1MWxW7rkltKmdnWRlgDhxBsR0EIcdIRR1gpdqcsit0K5qKERBQUFBQcBEWwFRQUFBwERbAti6Nm5VHsTlkUuxXMQolhKygoKDgIioetoKCg4CAogq2goKDgICiCnUyEEO2FEBeFEPFCiPJJvhslhLghhLgqhGhkKxsNIYQYK4QIFkKcVS8f29omfQghGqv36Q0hxEhb22MKQojbQohA9X62zKwcVkAIMU8I8VgIcUFjXQYhxE4hxHX1q68tbfwQUQQ7+VwA2gD7NVcKIYoBHYHiQGNgumoGZrtlspSytHrZYmtjdKHeh9OAJkAxoJN6XzsSddT72Z77NC9AddxqMhLYLaUsCOxWf1ZIQRTBTiZSystSyqtavmoJLJdSxkgpg4AbQMWUtS5VUhG4IaW8JaV8DSxHta8VLIiUcj8QmmR1S2Ch+v1CoFVK2qSgCLY1yQ7c0/h8X73OXhkshDivvhW251tdR9uvSZHADiHEKSFEP1sbYyJZpZQPANSvWWxszweHkq3PCIQQu4BsWr76Vkq5XtdmWtbZrA+lvt8A/AOMR2XfeGAS0CvlrDMJu9qvZlBNShkihMgC7BRCXFF7swoKBlEE2wiklPXN2Ow+kFPjcw7AyLm/LY+xv0EIMRvYZGVzkoNd7VdTkVKGqF8fCyHWogrxOIpgPxJC+EkpHwgh/IDHtjboQ0MJiViPDUBHIYS7ECIvUBA4bmObtKI++RJojepBqr1yAigohMgrhHBD9WB3g41tMgohRBohhE/Ce6Ah9r2vk7IB6K5+3x3QdXepYCUUDzuZCCFaA38BmYHNQoizUspGUsqLQogVwCUgFhgkpYyzpa16mCCEKI0qtHAb6G9Ta/QgpYwVQgwGtgPOwDwp5UUbm2UsWYG16vkmXYB/pZTbbGuSdoQQy4DaQCYhxH3ge+BXYIUQojdwF2hvOws/TJSh6QoKCgoOghISUVBQUHAQFMFWUFBQcBAUwVZQUFBwEBTBVlBQUHAQFMFWUFBQcBAUwVawe4QQnkKIfUIIZyFEbSGE2QN7hBDLhRAFLWmfgkJKoQi2giPQC1hjoX7s/wBfW6AeBYUURxFsBZshhKigTjjloR4FeFEIUUJL0S5oGVWn3v6MECKfOqf3QiHEDnXO6TZCiAnq3NPbhBCu6s0OAPWFEMqgMQWHQxFsBZshpTyBarjzj8AEYImU8p2h2urh5/mklLeTrK8KzABaSilvqVfnB5qiSgO6BNgrpSwJRKnXI6WMR5XqtpSVfpaCgtVQvAwFWzMOVX6QaGColu8zAc+TrCuKagbvhgnJlNRslVK+EUIEohq2njDsOxDIo1HuMeAPnEqu8QoKKYniYSvYmgyAN+ADeGj5PkrL+geoBL5MkvUxkOhFv5Fv8y7E865z4qGuV0HBoVAEW8HWzAJGA0uB35J+KaUMA5yFEJqi/RxViONnIURtM9osBDhKwigFhUQUwVawGUKIbkCslPJfVJngKggh6mopugOorrlCSvkIaA5ME0JUMqHNrEBUwswpCgqOhJKtT8HuEUKUAb6UUna1QF3DgJdSyrnJt0xBIWVRPGwFu0dKeQbYa6FZ55/zdiJZBQWHQvGwFRQUFBwExcNWUFBQcBAUwVZQUFBwEBTBVlBQUHAQFMFWUFBQcBAUwVZQUFBwEP4PkQ1xLgSHDyYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "this_y = 0 ## location that we want to take a slice through \n", " ## here setting constant y, so it's a west-east cross-section\n", " \n", "## bottom and top of the cross section in height; in km\n", "z_bot = 0\n", "z_top = 14 \n", "\n", "#for tt in range(0,num_times): ## loop over all the times\n", "for tt in range(6,7): ## just a specific time\n", "\n", " this_time = times[tt]\n", " this_time_minutes = this_time.dt.seconds.values/60.\n", " print(\"working on \"+str(this_time_minutes)+ \" minutes\")\n", " \n", " fig, ax = plt.subplots()\n", "\n", " ## color contoured vertical velocity\n", " cf = plt.contourf(y,heights.sel(zh=slice(z_bot,z_top)),\n", " www.sel(time=this_time,yh=this_y, method='nearest').sel(zh=slice(z_bot,z_top)),\n", " levels=np.arange(-30,33,3), cmap='bwr', extend='both')\n", " \n", " ## add pressure perturbations\n", " cf2 = plt.contour(y,heights.sel(zh=slice(z_bot,z_top)),\n", " ppert.sel(time=this_time,yh=this_y, method='nearest').sel(zh=slice(z_bot,z_top)),\n", " levels=np.arange(-250,275,25), colors='black', linewidths=1.0)\n", " \n", " plt.xlabel(\"x (km)\")\n", " plt.ylabel(\"height (km)\")\n", " \n", " cb = plt.colorbar(cf, ax=ax, shrink=0.9)\n", " cb.ax.tick_params(labelsize='x-small')\n", " \n", " plt.title(run_name+\" vertical velocity (m/s), ppert (Pa) y=\"+str(this_y)+\" time=\"+str(this_time_minutes)+ \" min\", fontsize='small')\n", " \n", " plt.savefig(run_name+\"/xsec_www_ppert_\"+\"{:04d}\".format(int(this_time_minutes))+\".png\", \n", " bbox_inches='tight',dpi=150, facecolor='white')\n", " \n", " ## uncomment this if you want to see the image in the notebook, comment it if you're doing a loop\n", " plt.show()\n", "\n", " plt.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }