From 0085a3a596d8c17d7e515b78b806c8d5dbf2ff69 Mon Sep 17 00:00:00 2001 From: Ellie Daw Date: Tue, 1 Oct 2024 10:51:13 -0500 Subject: [PATCH 1/3] Cocip grid example notebook --- docs/notebooks/CoCiPGrid.ipynb | 717 +++++++++++++++++++++++++++++++++ 1 file changed, 717 insertions(+) create mode 100644 docs/notebooks/CoCiPGrid.ipynb diff --git a/docs/notebooks/CoCiPGrid.ipynb b/docs/notebooks/CoCiPGrid.ipynb new file mode 100644 index 000000000..d92de6243 --- /dev/null +++ b/docs/notebooks/CoCiPGrid.ipynb @@ -0,0 +1,717 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from pycontrails.models.cocipgrid import CocipGrid\n", + "from pycontrails.datalib.ecmwf import ERA5\n", + "from pycontrails.core import MetDataset\n", + "import pandas as pd\n", + "from pycontrails.models.ps_model import PSGrid\n", + "from pycontrails.models.humidity_scaling import HistogramMatching" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up time and spatial bounds for the model run\n", + "time_bounds = (\"2022-03-01\", \"2022-03-02\")\n", + "lon_bounds = (-115, -95)\n", + "lat_bounds = (30, 50)\n", + "pressure_levels = [250, 300] # hPa" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [], + "source": [ + "# Download meteorological data\n", + "era5 = ERA5(time_bounds, pressure_levels=pressure_levels, variables=CocipGrid.met_variables)\n", + "met = era5.open_metdataset()\n", + "\n", + "era5_rad = ERA5(time_bounds, variables=CocipGrid.rad_variables)\n", + "rad = era5_rad.open_metdataset()" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/n1/l_d2y6h12xqcs35n7tpb6b_w0000gn/T/ipykernel_95121/1280914107.py:8: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.\n", + " \"time\": pd.date_range(time_bounds[0], time_bounds[1], freq=\"3H\"),\n" + ] + } + ], + "source": [ + "# Initialize CocipGrid model\n", + "ps_grid = PSGrid()\n", + "cocip_grid = CocipGrid(met, rad, aircraft_performance=ps_grid, humidity_scaling=HistogramMatching())\n", + "\n", + "# Define grid parameters\n", + "grid_params = {\n", + " \"level\": pressure_levels,\n", + " \"time\": pd.date_range(time_bounds[0], time_bounds[1], freq=\"3H\"),\n", + " \"longitude\": np.arange(lon_bounds[0], lon_bounds[1], 1.0),\n", + " \"latitude\": np.arange(lat_bounds[0], lat_bounds[1], 1.0),\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/ellie/devspace/breakthrough_analysis/venv/lib/python3.10/site-packages/pycontrails/models/cocipgrid/cocip_grid.py:2426: UserWarning: Met data 'rad' does not cover the source domain along the time axis. This causes some interpolated values to be nan, leading to meaningless results.\n", + " warnings.warn(\n", + "/Users/ellie/devspace/breakthrough_analysis/venv/lib/python3.10/site-packages/pycontrails/models/cocipgrid/cocip_grid.py:2568: UserWarning: End time of parameter 'met' (2022-03-02 00:00:00) is before model end time (2022-03-02 20:30:00). Include additional time at the end of 'met' or reduce 'max_age' parameter.\n", + " warnings.warn(\n", + "/Users/ellie/devspace/breakthrough_analysis/venv/lib/python3.10/site-packages/pycontrails/models/cocipgrid/cocip_grid.py:2568: UserWarning: End time of parameter 'rad' (2022-03-01 23:30:00) is before model end time (2022-03-02 20:30:00). Include additional time at the end of 'rad' or reduce 'max_age' parameter. Note: differencing reduces time coverage when providing accumulated radiative fluxes.\n", + " warnings.warn(\n", + "CocipGrid eval: 73%|███████▎ | 65/89 [00:04<00:01, 15.73it/s]\n" + ] + } + ], + "source": [ + "# Create a grid source\n", + "grid_source = MetDataset.from_coords(**grid_params)\n", + "\n", + "# Run CocipGrid model\n", + "result = cocip_grid.eval(source=grid_source)" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8QAAAK9CAYAAADi9qVKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACOdklEQVR4nOzdeZyNdeP/8fdhVsOM3VjG1mCIkYgGoSxjuUUqYTS2KGkjWcpOiEiliNta3CrFt7vbkmSojJ2Gui0jGtsQMmMsg5nz+6PfnNsx25kzZ876ej4e55Hzua7rc32u65o5p/d8PtfnMhiNRqMAAAAAAPAwhRzdAAAAAAAAHIFADAAAAADwSARiAAAAAIBHIhADAAAAADwSgRgAAAAA4JEIxAAAAAAAj0QgBgAAAAB4JAIxAAAAAMAjEYgBAAAAAB6JQAwALqxv376qWrWqo5shSWrVqpVatWqV63oxMTEyGAyKiYkp8DYh/5YuXSqDwaCTJ086uikAANgcgRgAbOT48eN6/vnnVb16dfn5+SkwMFDNmjXT+++/rxs3bji6edmKi4tTv379VK1aNfn5+alo0aJ64IEHNGLECP3+++8Oa1dGEMvutWPHDoe1zR6yO+7g4GBHNw0AALfh5egGAIA7+M9//qOnn35avr6+io6OVt26dXXr1i399NNPeuONN/Trr79qwYIFNt/vwoULlZ6enq/tBw8erNKlSysqKkphYWG6c+eODh06pOXLl2vOnDm6ceOGChcunGtd3333ndXtyMmkSZNUrVq1TOWhoaEFsj9n0rZtW0VHR5uV+fv727UNzz77rHr06CFfX1+77hcAAHsgEANAPp04cUI9evRQlSpV9MMPP6h8+fKmZUOGDFF8fLz+85//FMi+vb29rd52+/btGjx4sJo1a6Zvv/1WxYoVM1s+a9Ysvf3227nWc/36dRUpUkQ+Pj5WtyUnHTp0UKNGjQqk7ry4du2aAgIC7LrPmjVrqnfv3jav986dO0pPT7fomhUuXNiiP4gAAOCKGDINAPk0Y8YMpaSkaNGiRWZhOENoaKheffVV0/s7d+5o8uTJuu++++Tr66uqVavqzTffVGpqaqZt169fr5YtW6pYsWIKDAzUQw89pJUrV5qW33sP8cmTJ2UwGPTuu+/qvffeU5UqVeTv76+WLVvq0KFDZnVPnDhRBoNBK1asyBSGJcnPz0+TJ082C0OtWrVS3bp1tXfvXrVo0UJFihTRm2++aVp27z3Ep0+fVteuXRUQEKCyZctq6NChWR5nftx9zAsWLDCd14ceeki7d+/OtP7hw4f11FNPqWTJkvLz81OjRo30zTffmK2TMVx769atevHFF1W2bFlVqlTJtPyjjz5S9erV5e/vr8aNG+vHH380O/6UlBQFBASYXfe7z0nhwoU1bdq0fB/7hQsXNGDAAJUrV05+fn6qX7++li1blu35mTNnjun8/Pbbb6bz0b17d5UpU0b+/v6qVauW3nrrrUzn4u57iKtWrap//OMf+umnn9S4cWP5+fmpevXqWr58eaY2xsXFqWXLlvL391elSpU0ZcoULVmyhPuSAQBOgR5iAMinf//736pevbqaNm1q0frPPfecli1bpqeeekqvv/66du7cqWnTpum///2v1qxZY1pv6dKl6t+/v+6//36NHj1axYsX1/79+7Vhwwb16tUrx30sX75cV69e1ZAhQ3Tz5k29//77euyxx3Tw4EGVK1dO169f1w8//KBWrVqZBT1LXLp0SR06dFCPHj3Uu3dvlStXLsv1bty4odatWyshIUGvvPKKKlSooE8//VQ//PBDnvaXlJSkixcvmpUZDAaVKlXKrGzlypW6evWqnn/+eRkMBs2YMUPdunXT77//bupJ//XXX9WsWTNVrFhRo0aNUkBAgL744gt17dpVX331lZ544gmzOl988UWVKVNG48aN07Vr1yRJ8+bN00svvaRHHnlEQ4cO1cmTJ9W1a1eVKFHCdC6LFi2qJ554Qp9//rlmz55t9keFf/3rXzIajYqKisr12G/evJnp2IsVKyZfX1/duHFDrVq1Unx8vF566SVVq1ZNX375pfr27asrV65kCuNLlizRzZs3NWjQIPn6+qpkyZKKi4vTI488Im9vbw0aNEhVq1bV8ePH9e9//zvX0QHx8fF66qmnNGDAAPXp00eLFy9W37591bBhQ91///2SpDNnzujRRx+VwWDQ6NGjFRAQoH/+858MvwYAOA8jAMBqSUlJRknGLl26WLT+gQMHjJKMzz33nFn58OHDjZKMP/zwg9FoNBqvXLliLFasmLFJkybGGzdumK2bnp5u+nefPn2MVapUMb0/ceKEUZLR39/fePr0aVP5zp07jZKMQ4cONRqNRuMvv/xilGR87bXXMrXx0qVLxj///NP0Sk1NNS1r2bKlUZJx/vz5mbZr2bKlsWXLlqb3c+bMMUoyfvHFF6aya9euGUNDQ42SjFu2bMnhTBmNS5YsMUrK8uXr65vpmEuVKmW8fPmyqfz//u//jJKM//73v01lrVu3NtarV8948+ZNU1l6erqxadOmxho1amTad/PmzY137twxlaemphpLlSplfOihh4y3b982lS9dutQoyez4N27caJRkXL9+vdlxhYeHm62XneyOfcmSJUaj8X/n97PPPjNtc+vWLWNERISxaNGixuTkZLPzExgYaLxw4YLZPlq0aGEsVqyY8Y8//jArv/tnLONcnDhxwlRWpUoVoyTjtm3bTGUXLlww+vr6Gl9//XVT2csvv2w0GAzG/fv3m8ouXbpkLFmyZKY6AQBwBIZMA0A+JCcnS1KWQ46zsm7dOknSsGHDzMpff/11STLda7xp0yZdvXpVo0aNkp+fn9m6BoMh1/107dpVFStWNL1v3LixmjRpYtp/RruLFi2aadvq1aurTJkypte9w4l9fX3Vr1+/XNuwbt06lS9fXk899ZSprEiRIho0aFCu297to48+0qZNm8xe69evz7TeM888oxIlSpjeP/LII5Jkmin78uXL+uGHH9S9e3ddvXpVFy9e1MWLF3Xp0iVFRkbq2LFjOnPmjFmdAwcONOvd3bNnjy5duqSBAwfKy+t/g6yioqLM9i1Jbdq0UYUKFbRixQpT2aFDhxQXF2fxfcFdunTJdOyRkZGS/j6/wcHB6tmzp2l9b29vvfLKK0pJSdHWrVvN6nryySdVpkwZ0/s///xT27ZtU//+/VW5cmWzdS35GatTp47pHEtSmTJlVKtWLbOZyTds2KCIiAg98MADprKSJUta1DsOAIA9MGQaAPIhMDBQknT16lWL1v/jjz9UqFChTDMkBwcHq3jx4vrjjz8k/f0IJ0mqW7euVe2qUaNGprKaNWvqiy++kPS/AJ+SkpJpvf/7v//T7du39csvv2j48OGZllesWNGiyZj++OMPhYaGZgpXtWrVsugYMjRu3NiiSbXuDXUZAfWvv/6S9PcQX6PRqLFjx2rs2LFZ1nHhwgWzPyTcO7t1xvW59/p5eXlleh50oUKFFBUVpXnz5pkmHluxYoX8/Pz09NNP53o8klSpUiW1adMmy2V//PGHatSooUKFzP+2Xbt2bbO2ZncsGcHV2p+xe8+39Pc5zzjfGW2IiIjItJ4nzBAOAJbYtm2bZs6cqb179+rcuXNas2aNunbtmqc6Nm7cqPHjx+vXX3+Vn5+fWrRooVmzZmX6XkLWCMQAkA+BgYGqUKFCpgmrcmNJD1xBCg0NlZeXV5btbtmypSSZ9YDezd6P/bFUdjMhG41GSTI9nmr48OGmXtZ73RvU8nus0dHRmjlzptauXauePXtq5cqV+sc//qGgoKB81WsNW1+33M43ACB3165dU/369dW/f39169Ytz9ufOHFCXbp00bBhw7RixQolJSVp6NCh6tatm/bt21cALXY/BGIAyKd//OMfWrBggWJjY7PsDbtblSpVlJ6ermPHjpl68iTp/PnzunLliqpUqSJJuu+++yT9PcTWmt60Y8eOZSo7evSo6a/FAQEBatWqlbZu3aozZ86Y9YraSpUqVXTo0CEZjUazPwAcOXLE5vuyRPXq1SX9Paw4u17X3GRcn/j4eD366KOm8jt37ujkyZMKDw83W79u3bpq0KCBVqxYoUqVKikhIUEffvihlUeQuS1xcXFKT0836yU+fPiwWVuzk3E+8vrHnLy2MT4+PlN5VmUA4Ik6dOigDh06ZLs8NTVVb731lv71r3/pypUrqlu3rt555x3TUw327t2rtLQ0TZkyxfRdMHz4cHXp0kW3b9/O1+MZPQX3EANAPo0YMUIBAQF67rnndP78+UzLjx8/rvfff1+S1LFjR0nSnDlzzNaZPXu2JKlTp06SpHbt2qlYsWKaNm2abt68abauJT1wa9euNbsfdteuXdq5c6fZl+64ceOUlpam3r17Zzl0Or89fR07dtTZs2e1evVqU9n169e1YMGCfNVrrbJly6pVq1b65JNPdO7cuUzL//zzz1zraNSokUqVKqWFCxfqzp07pvIVK1aYDRW+27PPPqvvvvtOc+bMUalSpXL8H5+86NixoxITE/X555+byu7cuaMPP/xQRYsWNfX0Z6dMmTJq0aKFFi9erISEBLNlturljYyMVGxsrA4cOGAqu3z5stl91QCA7L300kuKjY3VqlWrFBcXp6efflrt27c3/eG7YcOGKlSokJYsWaK0tDQlJSXp008/VZs2bQjDFqKHGADy6b777tPKlSv1zDPPqHbt2oqOjlbdunV169Ytbd++3fQoHEmqX7+++vTpowULFujKlStq2bKldu3apWXLlqlr166mXsfAwEC99957eu655/TQQw+pV69eKlGihH755Rddv34907Nm7xUaGqrmzZtr8ODBSk1NNYWxESNGmNZ55JFHNHfuXL388suqUaOGoqKiFBYWplu3buno0aNasWKFfHx8FBwcbNV5GThwoObOnavo6Gjt3btX5cuX16effqoiRYrkqZ7169ebej3v1rRpU1Mvp6U++ugjNW/eXPXq1dPAgQNVvXp1nT9/XrGxsTp9+rR++eWXHLf38fHRhAkT9PLLL+uxxx5T9+7ddfLkSS1dulT33XdflkPhe/XqpREjRmjNmjUaPHiwzf4HZdCgQfrkk0/Ut29f7d27V1WrVtXq1av1888/a86cORZN9PbBBx+oefPmevDBBzVo0CBVq1ZNJ0+e1H/+8x+zEGutESNG6LPPPlPbtm318ssvmx67VLlyZV2+fNnhtw4AgDNLSEjQkiVLlJCQoAoVKkj6u/d3w4YNWrJkiaZOnapq1arpu+++U/fu3fX8888rLS1NERERpkk0kTsCMQDYwOOPP664uDjNnDlT//d//6d58+bJ19dX4eHhmjVrlgYOHGha95///KeqV6+upUuXas2aNQoODtbo0aM1fvx4szoHDBigsmXLavr06Zo8ebK8vb0VFhamoUOH5tqe6OhoFSpUSHPmzNGFCxfUuHFjzZ07V+XLlzdbb/DgwYqIiNB7772nL7/8UomJifL29tZ9992nPn36aPDgwabh23lVpEgRbd68WS+//LI+/PBDFSlSRFFRUerQoYPat29vcT3jxo3LsnzJkiV5DsR16tTRnj17NHHiRC1dulSXLl1S2bJl1aBBg2z3c6+XXnpJRqNRs2bN0vDhw1W/fn198803euWVVzLNCC5J5cqVU7t27bRu3To9++yzeWpvTvz9/RUTE6NRo0Zp2bJlSk5OVq1atbRkyRLTH2ByU79+fe3YsUNjx47VvHnzdPPmTVWpUkXdu3e3SRtDQkK0ZcsWvfLKK5o6darKlCmjIUOGKCAgINvzBQD428GDB5WWlqaaNWualaempqpUqVKSpMTERA0cOFB9+vRRz549dfXqVY0bN05PPfWUNm3axB8eLWAwMvsFALiNkydPqlq1apo5c2aWM0SjYKSnp6tMmTLq1q2bFi5cmGn5E088oYMHD3Lv7P/32muv6ZNPPlFKSkq2k3MBgKcxGAxms0x//vnnioqK0q+//prps7Jo0aIKDg7W2LFjtWHDBu3evdu07PTp0woJCVFsbKwefvhhex6CS6KHGACAPLh586Z8fX3N/uq+fPlyXb582TTJyd3OnTun//znP3rrrbfs2ErncePGDbMZri9duqRPP/1UzZs3JwwDQA4aNGigtLQ0Xbhwwey573e7fv16psfvZXy2ZjxdATkjEAMAkAc7duzQ0KFD9fTTT6tUqVLat2+fFi1apLp165o9X/jEiRP6+eef9c9//lPe3t56/vnnHdhqx4mIiFCrVq1Uu3ZtnT9/XosWLVJycnK2z4IGAE+SkpJiNnroxIkTOnDggEqWLKmaNWsqKipK0dHRmjVrlho0aKA///xTmzdvVnh4uDp16qROnTrpvffe06RJk0xDpt98801VqVJFDRo0cOCRuQ4CMQAAeVC1alWFhITogw8+0OXLl1WyZElFR0dr+vTp8vHxMa23detW9evXT5UrV9ayZcusnpzM1XXs2FGrV6/WggULZDAY9OCDD2rRokVq0aKFo5sGAA63Z88es8f4DRs2TJLUp08fLV26VEuWLNGUKVP0+uuv68yZMypdurQefvhh/eMf/5AkPfbYY1q5cqVmzJihGTNmqEiRIoqIiNCGDRts/vx5d8U9xAAAAAAAnTlzRiNHjtT69et1/fp1hYaGasmSJWrUqFGu2/78889q2bKl6tata5MnFdgLPcQAAAAA4OH++usvNWvWTI8++qjWr1+vMmXK6NixYypRokSu2165ckXR0dFq3bq1zp8/b4fW2g49xAAAAADg4UaNGqWff/5ZP/74Y5637dGjh2rUqKHChQtr7dq19BC7uvT0dJ09e1bFihXj2V0AAACAgxmNRl29elUVKlTINKuys7t586Zu3brlkH0bjcZMecbX11e+vr6Z1v3mm28UGRmpp59+Wlu3blXFihX14osvauDAgTnuY8mSJfr999/12WefacqUKTZtvz0QiLNw9uxZhYSEOLoZAAAAAO5y6tQpVapUydHNsNjNmzdVKTBQl27fdsj+ixYtqpSUFLOy8ePHa8KECZnW/f333zVv3jwNGzZMb775pnbv3q1XXnlFPj4+6tOnT5b1Hzt2TKNGjdKPP/4oLy/XjJau2eoCVqxYMUnSYz02y8unqINb43qq31/wH1I1Q53zuhyNT8l9JRfz+6+nHd0EAADg4e7cStEPq1qb/j/dVdy6dUuXbt/WNw0eUICdn71+LS1Nj+8/oFOnTikwMNBUnlXvsPT3KNlGjRpp6tSpkv5+DvKhQ4c0f/78LANxWlqaevXqpYkTJ6pmzZoFcxB2QCDOQsawAi+fovImEOfZqWNXJEmh4ZULbB9/nJHCajrPtTl89O8g7OsfmMuariM+LkGS+B0AAABOw1VvZwwoXFgBXvYNxBkCAwPNAnF2ypcvrzp16piV1a5dW1999VWW61+9elV79uzR/v379dJLL0n6O1QbjUZ5eXnpu+++02OPPZb/AyhgBGIUmPi4hAINxYePpjg8FGcEYXeTEYYBAADgGZo1a6YjR46YlR09elRVqlTJcv3AwEAdPHjQrOzjjz/WDz/8oNWrV6tatWoF1lZbIhCjQLlrKHbXICwRhgEAADzR0KFD1bRpU02dOlXdu3fXrl27tGDBAi1YsMC0zujRo3XmzBktX75chQoVUt26dc3qKFu2rPz8/DKVOzPXmqINLqmgA5a9wunhoymml7siDAMAAHimhx56SGvWrNG//vUv1a1bV5MnT9acOXMUFRVlWufcuXNKSHCv/1/kOcRZSE5OVlBQkNpF7+T+SRsqyJ7iDAXRW+zOAfhuhGEAAOCsbt9K0XfLmygpKcmi+2GdRUau2Nyood3vIb52J02t9+x1uXNmbwyZht1kBC5XGUJNEAYAAADcG0OmYXfOPoTa3YdF340wDAAAAE9GIIZDOGMo9qQgLBGGAQAAAAIxHMYZQrEnTJSVFcIwAAAAwD3EcDB7PJZJyjzZlqcF4AwEYQAAAOB/6CGGw8XHJditt9gTe4MzEIYBAAAAcwRiOA1nGELtrgjDAAAAQGYEYjgVgpvtcU4BAACArBGI4XQIcLbDuQQAAACyRyCGUyLI5Y897ssGAAAAXB2BGE6LUGcdzhkAAABgGQIxnB4Bz3KcKwAAAMByBGK4BIJe7jhHAAAAQN4QiOEyCHzZ49wAAAAAeefl6AYAeZER/ELDKzu4Jc6BIAwAAABYjx5iuCSCIOcAAAAAyC8CMVyWJwdCTz52AAAAwFYIxHBpnhgMPfGYAQAAgIJAIIbL86SA6EnHCgAAABQ0JtWCW3D3ybYIwgAAAIDt0UMMt+KOwdEdjwkAAABwBgRiuB13CpDudCwAAACAsyEQwy25Q5B0h2MAAAAAnBmBGG7LlQOlK7cdAAAAcBVMqgW35mqTbRGEAQAAAPuhhxgewRWCpjO2MTS8sukFAAAAuBt6iOExnDFwOrN7Q3BoeGXOIQAAANwKPcQAMsmuR5jeYgAAALgTAjEAM5YEXkIxAAAA3AGBGICkvPf+EooBAADg6gjEAKwOt4RiAAAAuDICMeDh8htqCcUAAABwVQRiwIPZKswy2RYAAABcEYEY8FAFEWAJxQAAAHAlBGLAwxR0by6hGAAAAK6CQAx4EHuFVUIxAAAAXAGBGPAQ9g6phGIAAAA4O6cJxNOnT5fBYNBrr70mSTp58qQMBkOWry+//DLbevr27Ztp/fbt29vpKADn5KhwymRbAAAAcGZejm6AJO3evVuffPKJwsPDTWUhISE6d+6c2XoLFizQzJkz1aFDhxzra9++vZYsWWJ67+vra9sGAy7EGQJpaHhlxcclOLoZAAAAgBmHB+KUlBRFRUVp4cKFmjJliqm8cOHCCg4ONlt3zZo16t69u4oWLZpjnb6+vpm2BTyNMwThuxGKAQAA4GwcPmR6yJAh6tSpk9q0aZPjenv37tWBAwc0YMCAXOuMiYlR2bJlVatWLQ0ePFiXLl3Kcf3U1FQlJyebvQBX5mxhOIOztgsAAACeyaGBeNWqVdq3b5+mTZuW67qLFi1S7dq11bRp0xzXa9++vZYvX67NmzfrnXfe0datW9WhQwelpaVlu820adMUFBRkeoWEhOT5WABn4eyh09nbBwAAAM/hsCHTp06d0quvvqpNmzbJz88vx3Vv3LihlStXauzYsbnW26NHD9O/69Wrp/DwcN13332KiYlR69ats9xm9OjRGjZsmOl9cnIyoRguyVXCZkY7GUINAAAAR3JYD/HevXt14cIFPfjgg/Ly8pKXl5e2bt2qDz74QF5eXmY9uqtXr9b169cVHR2d5/1Ur15dpUuXVnx8fLbr+Pr6KjAw0OwFuBpXCcN3c8U2AwAAwH04rIe4devWOnjwoFlZv379FBYWppEjR6pw4cKm8kWLFunxxx9XmTJl8ryf06dP69KlSypfvny+2ww4I1cPlUy2BQAAAEdxWA9xsWLFVLduXbNXQECASpUqpbp165rWi4+P17Zt2/Tcc89lWU9YWJjWrFkj6e8Zq9944w3t2LFDJ0+e1ObNm9WlSxeFhoYqMjLSLscF2JOrh+EM7nIcAAAAcC0On2U6N4sXL1alSpXUrl27LJcfOXJESUlJkv5+VFNcXJwef/xx1axZUwMGDFDDhg31448/8ixiuB13C5HudjwAAABwfgaj0Wh0dCOcTXJysoKCgtQueqe8fXJ+5jHgCO4eHhlCDQAA7nb7Voq+W95ESUlJLjXfT0au2NyooQK8Cue+gQ1du5Om1nv2utw5szen7yEG8D+h4ZXdPgxL7h/4AQAA4BwIxICL8LSQ6GnHCwAAAPsjEAMuwFPDoaceNwAAAOyDQAw4OU8PhZ4yTBwAAAD2RyAGnBhB8H84FwAAALA1AjHghOgVzRrnBAAAALZEIAacDKEvZ5wfAAAA2AqBGHAihD3LcJ4AAABgCwRiwEkQ8vKGYeUAAADILwIx4AQIdgAAAHCkCRMmyGAwmL3CwsJy3ObKlSsaMmSIypcvL19fX9WsWVPr1q2zU4ttw8vRDQAAAAAAON7999+v77//3vTeyyv7uHjr1i21bdtWZcuW1erVq1WxYkX98ccfKl68uB1aajsEYgAAAACAvLy8FBwcbNG6ixcv1uXLl7V9+3Z5e3tLkqpWrVqArSsYDJkGAAAAADeVnJxs9kpNTc123WPHjqlChQqqXr26oqKilJCQkO2633zzjSIiIjRkyBCVK1dOdevW1dSpU5WWllYQh1Fg6CEGAAAAgAIUEnGfivl623WfV1NvS3v2KiQkxKx8/PjxmjBhQqb1mzRpoqVLl6pWrVo6d+6cJk6cqEceeUSHDh1SsWLFMq3/+++/64cfflBUVJTWrVun+Ph4vfjii7p9+7bGjx9fUIdlcwRiAAAAAHBTp06dUmBgoOm9r69vlut16NDB9O/w8HA1adJEVapU0RdffKEBAwZkWj89PV1ly5bVggULVLhwYTVs2FBnzpzRzJkzCcQAAAAAAMcLDAw0C8SWKl68uGrWrKn4+Pgsl5cvX17e3t4qXLiwqax27dpKTEzUrVu35OPjY3Wb7Yl7iAEAAAAAZlJSUnT8+HGVL18+y+XNmjVTfHy80tPTTWVHjx5V+fLlXSYMSwRiAAAAAPB4w4cP19atW3Xy5Elt375dTzzxhAoXLqyePXtKkqKjozV69GjT+oMHD9bly5f16quv6ujRo/rPf/6jqVOnasiQIY46BKswZBoAAAAAPNzp06fVs2dPXbp0SWXKlFHz5s21Y8cOlSlTRpKUkJCgQoX+158aEhKijRs3aujQoQoPD1fFihX16quvauTIkY46BKsQiAEAAADAw61atSrH5TExMZnKIiIitGPHjgJqkX0wZBoAAAAA4JEIxAAAAAAAj0QgBgAAAAB4JAIxAAAAAMAjEYgBAAAAAB6JQAwAAAAA8EgEYgAAAACARyIQAwAAAAA8EoEYAAAAAOCRCMQAAAAAAI9EIAYAAAAAeCQCMQAAAADAIxGIAQAAAAAeiUAMAAAAAPBIBGIAAAAAgEciEAMAAAAAPBKBGAAAAADgkQjEAAAAAACPRCAGAAAAAHgkAjEAAAAAwCMRiAEAAAAAHsnL0Q0AAACA6wkNr+zoJris+LgERzcBwP9HDzEAAADyhDCcP5w/wHkQiAEnwF+KrcN5AwD7Cg2vTJizEc4j4BwIxICTiI9LIODlAecKAOyLAGd7/IEBcDwCMeBkCHq54xwBgH0R2goW5xdwHAIx4IQIfFmjFx0A7I+wZh+cZ8AxCMSAkyL4meN8AID9EdLsi/MN2B+BGHBihMC/cR4AwL64t9VxOO+AfRGIASfn6cOEPfnYAcARCGSOxx8kAPshEAMuwhODoSceMwA4EiHMuXA9gIJHIAZciKcERE/vFQcARyB8OSeuC1CwCMSAi3H3oOjuxwcAzojQ5dy4PkDBIRADLshdQ6O7HhcAOCvuVXUdXCugYBCIARflbsOK3elYAMAVEK5cE9cNsC0CMeDi3CFIusMxAIArIVS5Nq4fYDsEYsANuGqgdLdebgBwBYQp98B1BGyDQAy4CVcLlq7WXgBwB4Qo98L1BPKPQAy4EVcJma7STgBwF0zI5L64tkD+EIgBN+Psw5CduW0A4I4IS56B6wxYh0AMuClnC57OHtQBwB0RkjwL1xvIOwIx4MacJYA6SzsAwJMQjjwT1x3IGwIx4OYcHUYdvX8A8DTcUwquP2A5AjHgARw1XJkwDAD2RRBCBv4wAliGQAx4EHsGVMIwANgX4QdZ4ecCyBmBGPAwBR1UmTwLAOyP0IOc8PMBZI9ADHigggqsBGEAsD/CDizBzwmQNQIx4KFsHV4JwwBgX9wjirzi5wXIjEAMeDBbDW8mDAOAfRFsYC3+kAKYIxADyFegJQwDgH0RZmAL/BwBfyMQA5CU92DL5FkAYH+EGNgSP0+AEwXi6dOny2Aw6LXXXjOVtWrVSgaDwez1wgsv5FiP0WjUuHHjVL58efn7+6tNmzY6duxYAbcecA+WBlyCMADYH+EFBYGfK3g6pwjEu3fv1ieffKLw8PBMywYOHKhz586ZXjNmzMixrhkzZuiDDz7Q/PnztXPnTgUEBCgyMlI3b94sqOYDbiW3sEsYBgD74p5PFDR+vuDJHB6IU1JSFBUVpYULF6pEiRKZlhcpUkTBwcGmV2BgYLZ1GY1GzZkzR2PGjFGXLl0UHh6u5cuX6+zZs1q7dm2226Wmpio5OdnsBXiy7IZDE4YBwL4IKrAX/vACT+XwQDxkyBB16tRJbdq0yXL5ihUrVLp0adWtW1ejR4/W9evXs63rxIkTSkxMNKsrKChITZo0UWxsbLbbTZs2TUFBQaZXSEiI9QcEuJG7AzBhGAAAAO7Gy5E7X7Vqlfbt26fdu3dnubxXr16qUqWKKlSooLi4OI0cOVJHjhzR119/neX6iYmJkqRy5cqZlZcrV860LCujR4/WsGHDTO+Tk5MJxfB49/6VODS8MqEYAAAAbsVhgfjUqVN69dVXtWnTJvn5+WW5zqBBg0z/rlevnsqXL6/WrVvr+PHjuu+++2zWFl9fX/n6+tqsPsDVZTdkilAMAAAAd+KwIdN79+7VhQsX9OCDD8rLy0teXl7aunWrPvjgA3l5eSktLS3TNk2aNJEkxcfHZ1lncHCwJOn8+fNm5efPnzctA5Cz3O4f4v4iAAAAuAuHBeLWrVvr4MGDOnDggOnVqFEjRUVF6cCBAypcuHCmbQ4cOCBJKl++fJZ1VqtWTcHBwdq8ebOpLDk5WTt37lRERESBHAfgTiwNu0y8AQAA4N6yeixuVubMmaNatWrJ399fISEhGjp0qEs94cdhQ6aLFSumunXrmpUFBASoVKlSqlu3ro4fP66VK1eqY8eOKlWqlOLi4jR06FC1aNHC7PFMYWFhmjZtmp544gnTBZsyZYpq1KihatWqaezYsapQoYK6du1q5yMEXIs1AZch1AAAAO4np8fi3m3lypUaNWqUFi9erKZNm+ro0aPq27evDAaDZs+ebafW5o9DJ9XKiY+Pj77//nvNmTNH165dU0hIiJ588kmNGTPGbL0jR44oKSnJ9H7EiBG6du2aBg0apCtXrqh58+basGFDtvcpA54uvz29hGIAAAD3cfdjcadMmZLjutu3b1ezZs3Uq1cvSVLVqlXVs2dP7dy50x5NtQmnCsQxMTGmf4eEhGjr1q25bmM0Gs3eGwwGTZo0SZMmTbJ18wC3Y6thz4RiAAAA55ScnGz2PrcJhe9+LG5ugbhp06b67LPPtGvXLjVu3Fi///671q1bp2effdYmbbcHpwrEAOzH1vcAE4oBAACyFhgepsAi9h2xarj+93289z5Odvz48ZowYUKW2+T2WNx79erVSxcvXlTz5s1lNBp1584dvfDCC3rzzTfz1XZ7IhADHqigJsTKqJdgDAAA4BxOnTqlwMBA0/vseocteSzuvWJiYjR16lR9/PHHatKkieLj4/Xqq69q8uTJGjt2rE3aX9AIxIAHsdfM0PQWAwAAOIfAwECzQJydux+LmyEtLU3btm3T3LlzlZqamulJQGPHjtWzzz6r5557TpJUr14903xOb731lgoVcthDjSxGIAY8hL0fk0QoBgAAcB0Zj8W9W79+/RQWFqaRI0dm+Vjc69evZwq9GevdO9eTsyIQAx7AUc8MJhQDAAC4htweiytJ0dHRqlixoqZNmyZJ6ty5s2bPnq0GDRqYhkyPHTtWnTt3zjJAOyMCMeDmHBWG790/wRgAAMC1JSQkmPUIjxkzRgaDQWPGjNGZM2dUpkwZde7cWW+//bYDW5k3BGLAjTk6DN+N3mIAAADXcvdjcbN67+XlpfHjx2v8+PH2a5SNOf9dzgDyLDS8slOF4QzO2CYAAAB4LgIx4GacPXQ6e/sAAADgOQjEgBtxlbDpKu0EAACAeyMQA27C1UKmsw7rBgAAgOcgEANuwJWDpSu3HQAAAK6NQAy4MHfpZXWHYwAAAIDrIRADLsrdQqS7HQ8AAACcH4EYcEHuGh7d9bgAAADgnAjEgItx99DoLsPAAQAA4PwIxIAL8aSg6EnHCgAAAMcgEAMuwFN7TT3xmAEAAGA/BGLAyXl6KPT04wcAAEDBIRADToww+DfOAwAAAAoCgRhwUoRAc546bBwAAAAFh0AMOCGCX/Y4NwAAALAVAjHgROgFtQznCAAAALZAIAacBCEvbzhfAAAAyC8CMeAECHfW4bwBAAAgPwjEAAAAAACPRCAGnEB8XILi4xIc3QyXwzkD4O74nIM98fMGT0QgBpwIX0SW41wB8BR83sEe+DmDp/JydAMAmIuPS+De2BzwhQ3AE2V89vH9AFvjexWejh5iwAnx5ZQ1zgsAT8fnIGyJnyeAQAw4Lb6kzHE+AOBvfB7CFvg5Av5GIAacGJNt/Y1zAADm+FxEfvDzA/wPgRhwAZ78xeXJxw4AOeHzEXnFH9qBzAjEgIvwtC8wvrQBIHd8VsJS/JwAWSMQAy7EU77MPOU4AcBW+NxETvj5ALJHIAZcjLt/qbn78QFAQeHzE1nh5wLIGYEYcEHuOkTOHY8JAOyJz1HcjZ8HIHcEYsCFudMXnTsdCwA4Ep+ncNc/nAMFgUAMuDhX/8LjSxsAbI/PVs/FdQfyhkAMuAFX/fJz1XYDgKvgc9azcL2BvCMQA27C1b4EXa29AOCq+Lz1DFxnwDoEYsCNuMoQOVdoIwC4Ez533RvXF7AegRhwQ878xejMbQMAd8bnr/txlT+EA86MQAy4KWf7guRLGwAcj89i98F1BGyDQAy4MWf5snSWdgAA/sbnsmvj+gG2QyAG3JyjvzQdvX8AQNb4fHZNXDfAtgjEgAdw1BA5vrQBwLnxOe1auF6A7RGIAQ9iry9S7lEDANfB57Xz43sVKDgEYsDDFPQXKl/YAOB6CFzOi+sCFCwCMeCBCurLlS9tAHBtfI47F64HUPAIxICHsnVvAF/aAOAe+Dx3DlwHwD4IxICHs8UXLl/aAOBe+Fx3LM4/YD9ejm4AAMeLj0tQaHhlq7YDALinjM94a74fYB2+VwH7o4cYgKS8fwnzpQ0AnoHPe/vgPAOOQSAGYGLplzFf2gDgWfjcL1icX8BxCMQAzOQ22RZf2gDgmfj8LxicV8CxCMQAspTVFzRf2gDg2fgesB2e/Qw4BwIxgGxlfFHzpQ0AyMD3Qf454zkMDa9segGehEAMIEfO+KUNAHAcApP7ufeaco3hSQjEAAAAsAhByf1kd0251vAUBGIAAADkioDkfnK7pgyhhicgEAMAACBHhCL3k5dryvWHs7h9+7ZOnTqlI0eO6PLlyzapk0AMAACAbBGG3Iu1vb78HMBRrl69qnnz5qlly5YKDAxU1apVVbt2bZUpU0ZVqlTRwIEDtXv3bqvrJxADAAAgS4Qg95Lf68nPA+xt9uzZqlq1qpYsWaI2bdpo7dq1OnDggI4eParY2FiNHz9ed+7cUbt27dS+fXsdO3Ysz/sgEAMAACATwo97sdX15L5izzF9+nQZDAa99tprOa735ZdfKiwsTH5+fqpXr57WrVtnszbs3r1b27Zt065duzR27FhFRkaqXr16Cg0NVePGjdW/f38tWbJEiYmJ6tq1q3788cc874NADAAAADMEHvdSENeTnxH3tnv3bn3yyScKDw/Pcb3t27erZ8+eGjBggPbv36+uXbuqa9euOnTokE3a8a9//Uv3339/ruv5+vrqhRdeUP/+/fO8DwIxAAAATAg67qUgryc/K+4pJSVFUVFRWrhwoUqUKJHjuu+//77at2+vN954Q7Vr19bkyZP14IMPau7cuXZqbf55OboBAAAAcDzCjXux1/UMDa+s+LgEu+wL1klOTjZ77+vrK19f32zXHzJkiDp16qQ2bdpoypQpOdYdGxurYcOGmZVFRkZq7dq1Vrc3Ozdv3tSHH36oLVu26MKFC0pPTzdbvm/fPqvqJRADAAB4OMKwe7H39SQU5y6teh2lBRSx7z6vXZckhYSEmJWPHz9eEyZMyHKbVatWad++fRbP2pyYmKhy5cqZlZUrV06JiYl5b3AuBgwYoO+++05PPfWUGjduLIPBYJN6CcQAAAAejDDsXhx1PTP2SzB2PqdOnVJgYKDpfXa9w6dOndKrr76qTZs2yc/Pz17Ns9i3336rdevWqVmzZjat12nuIb53FrPLly/r5ZdfVq1ateTv76/KlSvrlVdeUVJSUo719O3bVwaDwezVvn17OxwBAACAayEMuxdnuJ7O0AaYCwwMNHtlF4j37t2rCxcu6MEHH5SXl5e8vLy0detWffDBB/Ly8lJaWlqmbYKDg3X+/HmzsvPnzys4ONjmx1GxYkUVK1bM5vU6RSDOahazs2fP6uzZs3r33Xd16NAhLV26VBs2bNCAAQNyra99+/Y6d+6c6fWvf/2rIJsPAADgcggu7sWZrqcztQWWa926tQ4ePKgDBw6YXo0aNVJUVJQOHDigwoULZ9omIiJCmzdvNivbtGmTIiIibN6+WbNmaeTIkfrjjz9sWq/Dh0zfPYvZ3Tdt161bV1999ZXp/X333ae3335bvXv31p07d+TllX3TfX198/RXidTUVKWmppre33vjOQAAgDshsLgPZ72W3FfseooVK6a6deualQUEBKhUqVKm8ujoaFWsWFHTpk2TJL366qtq2bKlZs2apU6dOmnVqlXas2ePFixYYPP2NWrUSDdv3lT16tVVpEgReXt7my2/fPmyVfU6PBDnZRazpKQkBQYG5hiGJSkmJkZly5ZViRIl9Nhjj2nKlCkqVapUtutPmzZNEydOtKr9AAAArsRZAxTyztmvJaHY/SQkJKhQof8NMm7atKlWrlypMWPG6M0331SNGjW0du3aTMHaFnr27KkzZ85o6tSpKleunHtMqpWXWcwuXryoyZMna9CgQTmu1759e3Xr1k3VqlXT8ePH9eabb6pDhw6KjY3NsptfkkaPHm02XXhycnKm2dgAAABcnbMHKFjOVa4lk225tpiYmBzfS9LTTz+tp59+usDbsn37dsXGxqp+/fo2rddhgTgvs5glJyerU6dOqlOnTrZThGfo0aOH6d/16tVTeHi47rvvPsXExKh169ZZbpPbs7gAAABcnasEKOTOFa8lvcXIr7CwMN24ccPm9TpsUi1LZzG7evWq2rdvr2LFimnNmjWZxornpnr16ipdurTi4+ML4jAAAACcnisGKGQWGl7Zpa+lK7cdjjd9+nS9/vrriomJ0aVLl5ScnGz2spbDeogzZjG7W79+/RQWFqaRI0eqcOHCSk5OVmRkpHx9ffXNN99Y9Tys06dP69KlSypfvrytmg4AAOAyCCHuwV2uIz3FsFbGo3TvHfVrNBplMBiyfCyUJRwWiHObxSw5OVnt2rXT9evX9dlnn5kl/zJlypjuBw4LC9O0adP0xBNPKCUlRRMnTtSTTz6p4OBgHT9+XCNGjFBoaKgiIyPtfowAAACO5C4hytO523UkFMMaP/zwg80m0rqbw2eZzs6+ffu0c+dOSVJoaKjZshMnTqhq1aqSpCNHjigpKUmSVLhwYcXFxWnZsmW6cuWKKlSooHbt2mny5MncIwwAADyKu4UoT+Wu15HJtmCpxYsX6/HHH1erVq0KpH6nCsR3z1rWqlUrGY3GXLe5ex1/f39t3LixIJoGAADgMtw1RHkaT7iO9BYjN5999plefPFFPfjgg+rSpYsef/xx1a5d22b1O2xSLQAAANieJ4Qod+fqk2fllScdK/Luhx9+0Llz5/Tiiy9q7969atKkiWrUqKHXX39d27ZtU3p6er7qJxADAAC4CYKF6/PUa+ipxw3LlChRQr1799YXX3yhixcv6sMPP9SNGzcUFRWlsmXLKjo6WqtXr9a1a9fyXDeBGAAAwA0QKFyfp19DTz9+WMbHx0ft27fXxx9/rFOnTmnDhg2qWrWqJk+erNmzZ+e5Pqe6hxgAAAB5R5BwfVzDvzHZFu61ZcsWNWvWTD4+Plkub9SokRo1aqRJkybp9u3bea6fQAwAAOCiCFHugeuYGZNtIUPr1q3l5+enhx9+WI8++qgeffRRPfzww/Lyyhxlvb2981w/Q6YBAABcECHK9Xna5Fl5xbmB9Pcjdz/66CNVrlxZixYtUosWLVS8eHFFRkZq+vTp2rlzZ74m1iIQAwAAuBiCguvjGlqG84QqVaqoX79+Wrp0qU6ePKn4+Hh98MEHKleunObNm6emTZuqZMmSVtfPkGkAAAAXQkBwfVzDvAkNr6z/7vnN0c2Ak6hevboKFy4sg8Egg8GgtWvX6tatW1bXRyAGAAAA4NSq31/J0U2AAyUkJCgmJkZbtmxRTEyMLl68qKZNm+qRRx7Rt99+qyZNmlhdN4EYAAAAAOCUqlevrr/++kvNmjVTixYt9Pzzz6tRo0ZZTqplDe4hBgAAcCHMvAvAk9y4cUOSVKhQIXl5ecnb21uFCxe2Wf0EYgAAABdDKAbgKc6dO6fY2Fh17NhRO3fuVKdOnVSiRAn94x//0Lvvvqvdu3czyzQAAICnIRQD8BRhYWF64YUX9PnnnysxMdEUkHft2qW2bdvma5ZpAjEAAICLIhQD8DTnz59XXFyc4uLi9Msvvyg5OVmpqalW18ekWgAAAC4sPi6Bx/gAcFsXLlxQTEyMaZbpo0ePytvbW40bN1aPHj306KOPKiIiwur6CcQAAAAujlAMwF0FBwfL29tbjRo10pNPPqlHH31UTZs2lb+/v03qJxADAAC4AUIxAHe0fv16NW/eXAEBAQVSP/cQAwAAuAnuKQbgbiIjIwssDEsEYgAAALdCKAbgLtq3b68dO3bkut7Vq1f1zjvv6KOPPsrzPhgyDQAA4GYyQrEnDqFm6DjgPp5++mk9+eSTCgoKUufOndWoUSNVqFBBfn5++uuvv/Tbb7/pp59+0rp169SpUyfNnDkzz/sgEAMAALgpTwuHGX8I8LTjBtzVgAED1Lt3b3355Zf6/PPPtWDBAiUlJUmSDAaD6tSpo8jISO3evVu1a9e2ah8EYgAAADfmCeGQYeKA+/L19VXv3r3Vu3dvSVJSUpJu3LihUqVKydvbO9/1cw8xAACAm3PXwBgfl5DtsbnrMQOeLigoyPQoJlsgEAMAAHgAdwuIlhyPux0zANsjEAMAAHgIdwiIOfUKA0BeEYgBAAA8iKuGSWuDsKseLwD7IBADAAB4GFcKifQIA0hLS9O2bdt05coVm9dNIAYAAPBArhAybdVGVzhWANkrXLiw2rVrp7/++svmdROIAQAAPJSzBkV6hQHcq27duvr9999tXi+BGAAAwIM5U/AsyCDsTMcJIO+mTJmi4cOH69tvv9W5c+eUnJxs9rKWlw3bCAAAABcUH5eg0PDKDm8DAGSnY8eOkqTHH39cBoPBVG40GmUwGJSWlmZVvQRiAAAAOCwU2zMIO0PwB2CdLVu2FEi9BGIAAABIsm9gpEcYQF60bNmyQOrlHmIAAACYFHRQdfSEWQRxwHX9+OOP6t27t5o2baozZ85Ikj799FP99NNPVtdJIAYAAIAZJrYC4Gy++uorRUZGyt/fX/v27VNqaqokKSkpSVOnTrW6XgIxAAAAMrFlT66je4Xv5UxtAWCZKVOmaP78+Vq4cKG8vb1N5c2aNdO+ffusrpdADAAAgGzlJzw6WxAG4LqOHDmiFi1aZCoPCgrSlStXrK6XQAwAAIAcWRNqnT0IO3v7AJgLDg5WfHx8pvKffvpJ1atXt7peAjEAAAByZWmApFcYQEEYOHCgXn31Ve3cuVMGg0Fnz57VihUrNHz4cA0ePNjqennsEgAAACyS02OZXDEE81xiwHWMGjVK6enpat26ta5fv64WLVrI19dXw4cP18svv2x1vQRiAAAAWOzeEOmKQRiA6zEYDHrrrbf0xhtvKD4+XikpKapTp46KFi2ar3oZMg0AAIA8yRgW7Q5h2B2OAfAE/fv319WrV+Xj46M6deqocePGKlq0qK5du6b+/ftbXS+BGAAAAB6NUAw4v2XLlunGjRuZym/cuKHly5dbXS9DpgEAAODxuJ8YcE7JyckyGo0yGo26evWq/Pz8TMvS0tK0bt06lS1b1ur6CcQAAAAAAKdUvHhxGQwGGQwG1axZM9Nyg8GgiRMnWl0/gRgAAAAQvcSAM9qyZYuMRqMee+wxffXVVypZsqRpmY+Pj6pUqaIKFSpYXT+BGAAAAPj/CMWAc2nZsqUk6cSJE6pcubIMBoNN62dSLQAAAOAuTLIFOJ8qVarop59+Uu/evdW0aVOdOXNGkvTpp5/qp59+srpeAjEAAABwD0Ix4Fy++uorRUZGyt/fX/v27VNqaqokKSkpSVOnTrW6XgIxAAAAkAVCMeA8pkyZovnz52vhwoXy9vY2lTdr1kz79u2zul4CMQAAAADAqR05ckQtWrTIVB4UFKQrV65YXS+BGAAAAMiGJ/USh9Us6ugmANkKDg5WfHx8pvKffvpJ1atXt7peAjEAAACQA08KxfBc8+bNU3h4uAIDAxUYGKiIiAitX78+2/UXLlyoRx55RCVKlFCJEiXUpk0b7dq1q8DaN3DgQL366qvauXOnDAaDzp49qxUrVmj48OEaPHiw1fXy2CUAAAAgF57yOKawmkV1+GiKo5sBB6hUqZKmT5+uGjVqyGg0atmyZerSpYv279+v+++/P9P6MTEx6tmzp5o2bSo/Pz+98847ateunX799VdVrFjR5u0bNWqU0tPT1bp1a12/fl0tWrSQr6+vhg8frpdfftnqeg1Go9Fow3a6heTkZAUFBald9E55+zB0BAAAAH9z51CcMWTaGQNx6o1kzR9ZUUlJSQoMDHR0cyyWkSvObF6lwIAi9t33teuq2LqHTp06ZXbOfH195evra1EdJUuW1MyZMzVgwIBc101LS1OJEiU0d+5cRUdHW93u3Ny6dUvx8fFKSUlRnTp1VLRo/vIaQ6YBAAAACzF8Gq4mJCREQUFBpte0adNy3SYtLU2rVq3StWvXFBERYdF+rl+/rtu3b6tkyZL5bXKOfHx8VKdOHTVu3DjfYVhiyDQAAACQJ+44fJoJtQrWnyVr6WaxYnbd51Xfq5KUZQ9xdg4ePKiIiAjdvHlTRYsW1Zo1a1SnTh2L9jdy5EhVqFBBbdq0yV/Ds3Hz5k19+OGH2rJliy5cuKD09HSz5dY+eolADAAAAABuKmOSLEvUqlVLBw4cUFJSklavXq0+ffpo69atuYbi6dOna9WqVYqJiZGfn58tmp3JgAED9N133+mpp55S48aNZTAYbFIvgRgAAADII3fsJQZ8fHwUGhoqSWrYsKF2796t999/X5988km227z77ruaPn26vv/+e4WHhxdY27799lutW7dOzZo1s2m93EMMAAAAWIH7ieHu0tPTlZqamu3yGTNmaPLkydqwYYMaNWpUoG2pWLGiihXAsHMCMQAAAGAldwzF3E/smUaPHq1t27bp5MmTOnjwoEaPHq2YmBhFRUVJkqKjozV69GjT+u+8847Gjh2rxYsXq2rVqkpMTFRiYqJSUgpmlvJZs2Zp5MiR+uOPP2xaL0OmAQAAgHxg+DTcwYULFxQdHa1z584pKChI4eHh2rhxo9q2bStJSkhIUKFC/+tPnTdvnm7duqWnnnrKrJ7x48drwoQJNm9fo0aNdPPmTVWvXl1FihSRt7e32fLLly9bVS+BGAAAAMgnQjFc3aJFi3JcHhMTY/b+5MmTBdeYLPTs2VNnzpzR1KlTVa5cOSbVAgAAAJB/DJGGK9i+fbtiY2NVv359m9bLPcQAAACADbjj/cSAswgLC9ONGzdsXq/Vgfj48eMaM2aMevbsqQsXLkiS1q9fr19//dVmjQMAAABcCaEYKBjTp0/X66+/rpiYGF26dEnJyclmL2tZFYi3bt2qevXqaefOnfr6669NM4n98ssvGj9+vNWNAQAAAFydO4RihlHD2bRv316xsbFq3bq1ypYtqxIlSqhEiRIqXry4SpQoYXW9Vt1DPGrUKE2ZMkXDhg0zexbUY489prlz51rdGAAAAMAdMMkWYFtbtmwpkHqtCsQHDx7UypUrM5WXLVtWFy9ezHejAAAAAFdHKAZsp2XLlgVSr1VDposXL65z585lKt+/f78qVqyY70YBAAAAAFDQrArEPXr00MiRI5WYmCiDwaD09HT9/PPPGj58uKKjo61qyPTp02UwGPTaa6+Zym7evKkhQ4aoVKlSKlq0qJ588kmdP38+x3qMRqPGjRun8uXLy9/fX23atNGxY8esahMAAACQH85+PzH3CsPTWRWIp06dqrCwMIWEhCglJUV16tRRixYt1LRpU40ZMybP9e3evVuffPKJwsPDzcqHDh2qf//73/ryyy+1detWnT17Vt26dcuxrhkzZuiDDz7Q/PnztXPnTgUEBCgyMlI3b97Mc7sAAACA/HL2UAx4MqsCsY+PjxYuXKjjx4/r22+/1WeffabDhw/r008/VeHChfNUV0pKiqKiorRw4UKz2cGSkpK0aNEizZ49W4899pgaNmyoJUuWaPv27dqxY0eWdRmNRs2ZM0djxoxRly5dFB4eruXLl+vs2bNau3atNYcKAAAA5BuhGLCe0WhUQkJCgXRyWv0cYkmqXLmyOnbsqO7du6tGjRpW1TFkyBB16tRJbdq0MSvfu3evbt++bVYeFhamypUrKzY2Nsu6Tpw4ocTERLNtgoKC1KRJk2y3kaTU1FSbPccKAAAAyIqrhWKGU8NZGI1GhYaG6tSpUzav2+JZpocNG2ZxpbNnz7ZovVWrVmnfvn3avXt3pmWJiYny8fFR8eLFzcrLlSunxMTELOvLKC9XrpzF20jStGnTNHHiRIvaDAAAAACwn0KFCqlGjRq6dOmS1R2x2bE4EO/fv9/s/b59+3Tnzh3VqlVLknT06FEVLlxYDRs2tKi+U6dO6dVXX9WmTZvk5+eXhybb3ujRo80Cf3JyskJCQhzYIgAAAABAhunTp+uNN97QvHnzVLduXZvVa3EgvvtByLNnz1axYsW0bNky032/f/31l/r166dHHnnEovr27t2rCxcu6MEHHzSVpaWladu2bZo7d642btyoW7du6cqVK2a9xOfPn1dwcHCWdWaUnz9/XuXLlzfb5oEHHsi2Lb6+vvL19bWo3QAAAIA7YEg0XEl0dLSuX7+u+vXry8fHR/7+/mbLL1++bFW9Fgfiu82aNUvfffed2SRYJUqU0JQpU9SuXTu9/vrrudbRunVrHTx40KysX79+CgsL08iRIxUSEiJvb29t3rxZTz75pCTpyJEjSkhIUERERJZ1VqtWTcHBwdq8ebMpACcnJ2vnzp0aPHiwNYcKAAAAAHCwOXPmFEi9VgXi5ORk/fnnn5nK//zzT129etWiOooVK5apqzsgIEClSpUylQ8YMEDDhg1TyZIlFRgYqJdfflkRERF6+OGHTduEhYVp2rRpeuKJJ0zPMZ4yZYpq1KihatWqaezYsapQoYK6du1qzaECAAAAABysT58+BVKvVYH4iSeeUL9+/TRr1iw1btxYkrRz50698cYbuT4nOC/ee+89FSpUSE8++aRSU1MVGRmpjz/+2GydI0eOKCkpyfR+xIgRunbtmgYNGqQrV66oefPm2rBhg8PvUwYAAABcTVjNojp8NMXRzQAkScePH9eSJUt0/Phxvf/++ypbtqzWr1+vypUr6/7777eqToPRaDTmdaPr169r+PDhWrx4sW7fvi1J8vLy0oABAzRz5kwFBARY1RhnkZycrKCgILWL3ilvH+6tAAAAgG2Ehld2dBNMLL2H2BkCceqNZM0fWVFJSUkKDAx0dHMslpErDuzfr2LFitl131evXtUDDRq43DnLztatW9WhQwc1a9ZM27Zt03//+19Vr15d06dP1549e7R69Wqr6rXqOcRFihTRxx9/rEuXLmn//v3av3+/Ll++rI8//tjlwzAAAAAAwLmMGjVKU6ZM0aZNm+Tj42Mqf+yxx7Rjxw6r67VqyHSGgIAAhYeH56cKAAAAAABydPDgQa1cuTJTedmyZXXx4kWr67UqED/66KMyGAzZLv/hhx+sbhAAAACAgsUjl+BqihcvrnPnzqlatWpm5fv371fFihWtrteqQHzvM31v376tAwcO6NChQwU2+xcAAAAAwDP16NFDI0eO1JdffimDwaD09HT9/PPPGj58uKKjo62u16pA/N5772VZPmHCBKWkOP6mewAAAACA+5g6daqGDBmikJAQpaWlqU6dOkpLS1OvXr00ZswYq+u1alKt7PTu3VuLFy+2ZZUAAAAAHIjh1XAGPj4+WrhwoY4fP65vv/1Wn332mQ4fPqxPP/1UhQsXtrrefE2qda/Y2Fie9wsAAAAAKBCVK1dWSEiIJOU4r5WlrArE3bp1M3tvNBp17tw57dmzR2PHjs13owAAAAAAuNuiRYv03nvv6dixY5KkGjVq6LXXXtNzzz1ndZ1WBeLAwECzNF6oUCHVqlVLkyZNUrt27axuDAAAAICCxRBouKJx48Zp9uzZevnllxURESHp7xHKQ4cOVUJCgiZNmmRVvVYF4qVLl1q1MwAAAAAA8mrevHlauHChevbsaSp7/PHHFR4erpdfftnqQGzVpFrVq1fXpUuXMpVfuXJF1atXt6ohAAAAAABk5fbt22rUqFGm8oYNG+rOnTtW12tVID558qTS0tIylaempurMmTNWNwYAAACA82GYNRzt2Wef1bx58zKVL1iwQFFRUVbXm6ch0998843p3xs3blRQUJDpfVpamjZv3qyqVata3RgAAAAAALKyaNEifffdd3r44YclSTt37lRCQoKio6M1bNgw03qzZ8+2uM48BeKuXbtK+nt66z59+pgt8/b2VtWqVTVr1qy8VAkAAAAAQI4OHTqkBx98UJJ0/PhxSVLp0qVVunRpHTp0yLReXh/FlKdAnJ6eLkmqVq2adu/erdKlS+dpZwAAAAAA5NWWLVsKpF6rZpk+ceKErdsBAAAAoIBxLzBgzuJA/MEHH2jQoEHy8/PTBx98kOO6r7zySr4bBgAAAABAQbI4EL/33nuKioqSn5+f3nvvvWzXMxgMBGIAAAAAgNOzOBDfPUyaIdMAAACAZwmrWVSHj6Y4uhmATVn1HOJJkybp+vXrmcpv3LihSZMm5btRAAAAAAAUNKsm1Zo4caJeeOEFFSlSxKz8+vXrmjhxosaNG2eTxgEAAAAA8M0332RZbjAY5Ofnp9DQUFWrVi3P9VoViI1GY5bPd/rll19UsmRJa6oEAAAAUICYYRqurGvXrjIYDDIajWblGWUGg0HNmzfX2rVrVaJECYvrzdOQ6RIlSqhkyZIyGAyqWbOmSpYsaXoFBQWpbdu26t69e16qBAAAAAAgR5s2bdJDDz2kTZs2KSkpSUlJSdq0aZOaNGmib7/9Vtu2bdOlS5c0fPjwPNWbpx7iOXPmyGg0qn///po4caKCgoJMy3x8fFS1alVFRETkqQEAAAAAAOTk1Vdf1YIFC9S0aVNTWevWreXn56dBgwbp119/1Zw5c9S/f/881ZunQNynTx9JUrVq1dS0aVN5e3vnaWcAAAAAAOTV8ePHFRgYmKk8MDBQv//+uySpRo0aunjxYp7qtWqW6ZYtW5rC8M2bN5WcnGz2AgAAAGAuNLyyo5sAuKyGDRvqjTfe0J9//mkq+/PPPzVixAg99NBDkqRjx44pJCQkT/VaFYivX7+ul156SWXLllVAQIBKlChh9gIAAADwP84QhnmGMFzZokWLdOLECVWqVEmhoaEKDQ1VpUqVdPLkSf3zn/+UJKWkpGjMmDF5qteqWabfeOMNbdmyRfPmzdOzzz6rjz76SGfOnNEnn3yi6dOnW1MlAAAA4JacIQzbAoEajlSrVi399ttv+u6773T06FFTWdu2bVWo0N/9vF27ds1zvVYF4n//+99avny5WrVqpX79+umRRx5RaGioqlSpohUrVigqKsqaagEAAAC34mxh+PDRFKsev0QYhjMoVKiQ2rdvr/bt29usTqsC8eXLl1W9enVJf9/EfPnyZUlS8+bNNXjwYJs1DgAAAIBtWRuKAUfbvHmzNm/erAsXLig9Pd1s2eLFi62q06p7iKtXr64TJ05IksLCwvTFF19I+rvn+O5HMQEAAACeytl6h61F7zCcwcSJE9WuXTtt3rxZFy9e1F9//WX2spZVPcT9+vXTL7/8opYtW2rUqFHq3Lmz5s6dq9u3b2v27NlWNwYAAABwB84ehi3tJSYMw1nMnz9fS5cu1bPPPmvTeq0KxEOHDjX9u02bNjp8+LD27t2r0qVL67PPPrNZ4wAAAABX4+xhOANDp+FKbt26paZNm9q8XquGTN+rSpUq6tatm4KCgrRo0SJbVAkAAAC4HFcJw5agdxjO5LnnntPKlSttXq9VPcQAAAAAzLliGKaXGK7i5s2bWrBggb7//nuFh4fL29vbbLm1t+4SiAEAAIB8csUwnCGrUEzvMJxNXFycHnjgAUnSoUOHzJYZDAar6yUQAwAAAPngymE4K4RhOKMtW7YUSL15CsTdunXLcfmVK1fy0xYAAADApbhLGGboNDxVngJxbs8YDgoKUnR0dL4aBAAAALgCdwnDGegZhrPp1q2bli5dqsDAwFw7Z7/++mur9pGnQLxkyRKrdgIAAAC4E3cLw8C0adP09ddf6/Dhw/L391fTpk31zjvvqFatWhZtv2rVKvXs2VNdunTR2rVrbdKmoKAg0/3BuXXOWot7iAEAAIA8IAzDHW3dulVDhgzRQw89pDt37ujNN99Uu3bt9NtvvykgICDHbU+ePKnhw4frkUcesWmb7u6QLajOWQIxAAAAYCHCMNzVhg0bzN4vXbpUZcuW1d69e9WiRYtst0tLS1NUVJQmTpyoH3/8scDmlTpx4oTu3LmjGjVqmJUfO3ZM3t7eqlq1qlX1FrJB2wAAAAAATig5OdnslZqaatF2SUlJkqSSJUvmuN6kSZNUtmxZDRgwIN9tzUnfvn21ffv2TOU7d+5U3759ra6XHmIAAADAAvQOw1p/3KiggMKBdt3ntRvJkqSQkBCz8vHjx2vChAk5bpuenq7XXntNzZo1U926dbNd76efftKiRYt04MCB/DY3V/v371ezZs0ylT/88MN66aWXrK6XQAwAAADkgjAMV3Xq1CkFBv4vjPv6+ua6zZAhQ3To0CH99NNP2a5z9epVPfvss1q4cKFKly5tk7bmxGAw6OrVq5nKk5KSlJaWZnW9BGIAAAAgB4RhuLLAwECzQJybl156Sd9++622bdumSpUqZbve8ePHdfLkSXXu3NlUlp6eLkny8vLSkSNHdN9991nf8Hu0aNFC06ZN07/+9S8VLlxY0t/3L0+bNk3Nmze3ul4CMQAAAJANwjA8hdFo1Msvv6w1a9YoJiZG1apVy3H9sLAwHTx40KxszJgxunr1qt5///1MQ7Xza/r06WrZsqVq1aplms36xx9/VHJysn744Qer6yUQAwAAAFkgDMOTDBkyRCtXrtT//d//qVixYkpMTJT09/N//f39JUnR0dGqWLGipk2bJj8/v0z3FxcvXlyScrzv2Fr333+/4uLiNHfuXP3yyy/y9/dXdHS0XnrppVwn/soJgRgAAAC4B2EYnmbevHmSpFatWpmVL1myxDSLc0JCggoVsv+Dim7fvq327dtr/vz5mjp1qk3rJhADAAAAdyEMwxMZjcZc14mJiclx+dKlS23TmHt4e3srLi6uQOrmOcQAAADA/0cYBpxT7969tWjRIpvXSw8xAAAAIMIw4Mzu3LmjxYsX6/vvv1fDhg0VEBBgtnz27NlW1UsgBgAAgMcjDAPO7dChQ3rwwQclSUePHjVbZjAYrK6XQAwAAACPRhgGnN+WLVsKpF7uIQYAAIDHIgwDruf06dM6ffq0TeoiEAMAAAAAnFp6eromTZqkoKAgValSRVWqVFHx4sU1efJkpaenW10vQ6YBAADgkegdBlzHW2+9pUWLFmn69Olq1qyZJOmnn37ShAkTdPPmTb399ttW1UsgBgAAgMchDAOuZdmyZfrnP/+pxx9/3FQWHh6uihUr6sUXX7Q6EDNkGgAAAADg1C5fvqywsLBM5WFhYbp8+bLV9RKIAQAAAABOrX79+po7d26m8rlz56p+/fpW18uQaQAAAHgUhksDrmfGjBnq1KmTvv/+e0VEREiSYmNjderUKa1bt87qeukhBgAAAAA4tZYtW+ro0aN64okndOXKFV25ckXdunXTkSNH9Mgjj1hdLz3EAAAAAACn9Pvvv6tatWoyGAyqUKGC1ZNnZYceYgAAAACAU6pRo4b+/PNP0/tnnnlG58+ft1n9BGIAAAAAgFMyGo1m79etW6dr167ZrH4CMQAAADwGE2oBuBuBGAAAAADglAwGgwwGQ6YyW2FSLQAAAACAUzIajerbt698fX0lSTdv3tQLL7yggIAAs/W+/vprq+onEAMAAAAAnFKfPn3M3vfu3dum9ROIAQAAAABOacmSJQVav0PvIZ43b57Cw8MVGBiowMBARUREaP369ZKkkydPmsaL3/v68ssvs62zb9++mdZv3769vQ4JAAAATooJtQDcy6E9xJUqVdL06dNVo0YNGY1GLVu2TF26dNH+/fsVFhamc+fOma2/YMECzZw5Ux06dMix3vbt25v9JSFjvDkAAAAAABkcGog7d+5s9v7tt9/WvHnztGPHDt1///0KDg42W75mzRp1795dRYsWzbFeX1/fTNsCAAAAAHA3p3nsUlpamlatWqVr164pIiIi0/K9e/fqwIEDGjBgQK51xcTEqGzZsqpVq5YGDx6sS5cu5bh+amqqkpOTzV4AAAAAAPfm8EB88OBBFS1aVL6+vnrhhRe0Zs0a1alTJ9N6ixYtUu3atdW0adMc62vfvr2WL1+uzZs365133tHWrVvVoUMHpaWlZbvNtGnTFBQUZHqFhITk+7gAAADgPLh/GHBt165dK5B6HR6Ia9WqpQMHDmjnzp0aPHiw+vTpo99++81snRs3bmjlypUW9Q736NFDjz/+uOrVq6euXbvq22+/1e7duxUTE5PtNqNHj1ZSUpLpderUqfweFgAAAADARsqVK6f+/fvrp59+smm9Dg/EPj4+Cg0NVcOGDTVt2jTVr19f77//vtk6q1ev1vXr1xUdHZ3n+qtXr67SpUsrPj4+23V8fX1NM11nvAAAAAAAzuGzzz7T5cuX9dhjj6lmzZqaPn26zp49m+96HR6I75Wenq7U1FSzskWLFunxxx9XmTJl8lzf6dOndenSJZUvX95WTQQAAAAA2FHXrl21du1anTlzRi+88IJWrlypKlWq6B//+Ie+/vpr3blzx6p6HRqIR48erW3btunkyZM6ePCgRo8erZiYGEVFRZnWiY+P17Zt2/Tcc89lWUdYWJjWrFkjSUpJSdEbb7yhHTt26OTJk9q8ebO6dOmi0NBQRUZG2uWYAAAAAAAFo0yZMho2bJji4uI0e/Zsff/993rqqadUoUIFjRs3TtevX89TfQ597NKFCxcUHR2tc+fOKSgoSOHh4dq4caPatm1rWmfx4sWqVKmS2rVrl2UdR44cUVJSkiSpcOHCiouL07Jly3TlyhVVqFBB7dq10+TJk3kWMQAAgIdiQi3AfZw/f17Lli3T0qVL9ccff+ipp57SgAEDdPr0ab3zzjvasWOHvvvuO4vrc2ggXrRoUa7rTJ06VVOnTs12udFoNP3b399fGzdutEnbAAAAAADO4euvv9aSJUu0ceNG1alTRy+++KJ69+6t4sWLm9Zp2rSpateunad6HRqIAQAAAADITb9+/dSjRw/9/PPPeuihh7Jcp0KFCnrrrbfyVC+BGAAAAADg1M6dO6ciRYrkuI6/v7/Gjx+fp3oJxAAAAAAAp3bnzh0lJydnKjcYDPL19ZWPj49V9RKIAQAA4LaYUAtwD8WLF5fBYMh2eaVKldS3b1+NHz9ehQpZ/jAlAjEAAAAAwKktXbpUb731lvr27avGjRtLknbt2qVly5ZpzJgx+vPPP/Xuu+/K19dXb775psX1EogBAAAAAE5t2bJlmjVrlrp3724q69y5s+rVq6dPPvlEmzdvVuXKlfX222/nKRBb3pcMAAAAAIADbN++XQ0aNMhU3qBBA8XGxkqSmjdvroSEhDzVSyAGAACAW+L+YcB9hISEaNGiRZnKFy1apJCQEEnSpUuXVKJEiTzVy5BpAAAAAIBTe/fdd/X0009r/fr1pucQ79mzR4cPH9bq1aslSbt379YzzzyTp3oJxAAAAAAUH5dArzqc1uOPP64jR47ok08+0ZEjRyRJHTp00Nq1a1W1alVJ0uDBg/NcL4EYAAAAbodgB7iP27dvq3379po/f76mTZtm07q5hxgAAAAA4LS8vb0VFxdXIHUTiAEAAOBW6B0G3E/v3r2znFQrvxgyDQAAALdBGAbc0507d7R48WJ9//33atiwoQICAsyWz54926p6CcQAAABwC4RhwH0dOnRIDz74oCTp6NGjZssMBoPV9RKIAQAA4NIIwrbBeYQz27JlS4HUyz3EAAAAcFmEOMCzxMfHa+PGjbpx44YkyWg05qs+AjEAAABcEmEY8ByXLl1S69atVbNmTXXs2FHnzp2TJA0YMECvv/661fUSiAEAAOBSQsMrE4YBDzN06FB5e3srISFBRYoUMZU/88wz2rBhg9X1cg8xAAAAXAZBGPBM3333nTZu3KhKlSqZldeoUUN//PGH1fUSiAEAAOD0CMKAZ7t27ZpZz3CGy5cvy9fX1+p6GTINAAAAp0YYBvDII49o+fLlpvcGg0Hp6emaMWOGHn30UavrpYcYAAAAToswDECSZsyYodatW2vPnj26deuWRowYoV9//VWXL1/Wzz//bHW9BGIAAAA4HYIwgLvVrVtXR48e1dy5c1WsWDGlpKSoW7duGjJkiMqXL291vQRiAAAAOBXCMICsBAUF6a233rJpnQRiAAAAOAWCMLJTM7Soo5sAJ3DlyhXt2rVLFy5cUHp6utmy6Ohoq+okEAMAAMDhCMPITljNorpxLdnRzXB727Zt08yZM7V3716dO3dOa9asUdeuXXPcJjU1VZMmTdJnn32mxMRElS9fXuPGjVP//v1t3r5///vfioqKUkpKigIDA2UwGEzLDAYDgRgAAACuhyCMnITVpGfYXq5du6b69eurf//+6tatm0XbdO/eXefPn9eiRYsUGhqqc+fOZeq5tZXXX39d/fv319SpU7N8/JK1CMQAAABwCMIwckIYtq8OHTqoQ4cOFq+/YcMGbd26Vb///rtKliwpSapatWoBtU46c+aMXnnlFZuGYYnnEAMAAMABCMOAfSQnJ5u9UlNTbVLvN998o0aNGmnGjBmqWLGiatasqeHDh+vGjRs2qf9ekZGR2rNnj83rpYcYAAAAdkMQhiXcrXf4+Dk/+Qf42XWfN67dkiSFhISYlY8fP14TJkzId/2///67fvrpJ/n5+WnNmjW6ePGiXnzxRV26dElLlizJd/336tSpk9544w399ttvqlevnry9vc2WP/7441bVSyAGAACAXRCGYQl3C8OOdurUKQUGBpre+/r62qTe9PR0GQwGrVixQkFBQZKk2bNn66mnntLHH38sf39/m+wnw8CBAyVJkyZNyrTMYDAoLS3NqnoJxAAAAChQBGFYijBse4GBgWaB2FbKly+vihUrmsKwJNWuXVtGo1GnT59WjRo1bLq/gpqsi3uIAQAAUGAIw7AUYdi1NGvWTGfPnlVKSoqp7OjRoypUqJAqVarkwJblDYEYAAAABYIwDEsRhh0vJSVFBw4c0IEDByRJJ06c0IEDB5SQkCBJGj16tNmzfnv16qVSpUqpX79++u2337Rt2za98cYb6t+/v02HS3fs2FFJSUmm99OnT9eVK1dM7y9duqQ6depYXT+BGAAAADYVGl6ZMAyLEYadw549e9SgQQM1aNBAkjRs2DA1aNBA48aNkySdO3fOFI4lqWjRotq0aZOuXLmiRo0aKSoqSp07d9YHH3xg03Zt3LjRbGbsqVOn6vLly6b3d+7c0ZEjR6yun3uIAQAAYDMEYeQFYdh5tGrVSkajMdvlS5cuzVQWFhamTZs2FWCrlKlNObXRGgRiAAAA5BtBGHlFGIYzYMg0AAAA8oUwjLwiDMNSBoNBBoMhU5mt0EMMAAAAqxGGARQko9Govn37mp6ffPPmTb3wwgsKCAiQJLP7i61BIAYAAECeEYRhLXqHkRd9+vQxe9+7d+9M69w9+3VeEYgBAAAA2AVhGHm1ZMmSAq2fe4gBAACQZ/FxCbmvBJdij2t6+GiKDh9NKfD9AJYiEAMAAMAqhGL3kXEt4+MS7BaMAWdAIAYAAIDV7BWgUHCyun6EYngKAjEAAADyjVDsmnK6boRieAICMQAAAGyCUOw6LO3ZJxTD3RGIAQAAYDMMoXZ+eb0+TLYFd0YgBgAAgM0Rip2TtdeFybbgrgjEAAAAKBCEYudii+tBKIa7IRADAACgwDCE2jnY8hoQiuFOCMQAAAAocIRixyioP0gQiuEuCMQAAACwC0KxfRX0+WayLbgDAjEAAADshiHU9mGvc2yv63k0nlCMgkEgBgAAgN0RiguOI84t1xOuikAMAAAAhyBE2Zaje9+5nnBFBGIAAAA4jKNDnLtwlnPoLO0ALEUgBgAAgMMRpKznbOfO2doD5IRADAAAAKdAkMo7Zz1n9PzDVRCIAQAA4DQIUpZzhfPkCm2EZyMQAwAAwOkQpLLnan80cKW2wvMQiAEAAOCUCFKZueo5cdV2w/0RiAEAAOC0XK03tCC5+nlw9fbDPRGIAQAA4PQ8PUy5y/HzBw44GwIxAAAAXIKnBil3PG53PCa4JgIxAAAAXIYn9TC6+7G687HBdRCIAQAA4HLcPUy5+/Fl8JTjhPMiEAMAAMAluWuYctfjyo6794TDuRGIAQAA4LLcLUy507HklScfOxyHQAwAAACX5w5hyh2OIb84B7A3hwbiefPmKTw8XIGBgQoMDFRERITWr19vWt6qVSsZDAaz1wsvvJBjnUajUePGjVP58uXl7++vNm3a6NixYwV9KAAAAHAwVw1T7tbLnV+cC9iTlyN3XqlSJU2fPl01atSQ0WjUsmXL1KVLF+3fv1/333+/JGngwIGaNGmSaZsiRYrkWOeMGTP0wQcfaNmyZapWrZrGjh2ryMhI/fbbb/Lz8yvQ4wEAAIBj5RamQsMr26klliH8ZS0+LsHprhXck0MDcefOnc3ev/3225o3b5527NhhCsRFihRRcHCwRfUZjUbNmTNHY8aMUZcuXSRJy5cvV7ly5bR27Vr16NHDtgcAAAAAl5LXAFqQoYwwnLOM80MwRkFyaCC+W1pamr788ktdu3ZNERERpvIVK1bos88+U3BwsDp37qyxY8dm20t84sQJJSYmqk2bNqayoKAgNWnSRLGxsdkG4tTUVKWmppreJycn2+ioAAAA4MqsCa2WBjhnDXrOFtTj4xIUUqO4o5sBN+XwQHzw4EFFRETo5s2bKlq0qNasWaM6depIknr16qUqVaqoQoUKiouL08iRI3XkyBF9/fXXWdaVmJgoSSpXrpxZebly5UzLsjJt2jRNnDjRRkcEAAAAT5ZdoHTWAHyv0PDKTheKf//1tKObADfl8EBcq1YtHThwQElJSVq9erX69OmjrVu3qk6dOho0aJBpvXr16ql8+fJq3bq1jh8/rvvuu89mbRg9erSGDRtmep+cnKyQkBCb1Q8AAADP5SpB+G7OGIqBguDwxy75+PgoNDRUDRs21LRp01S/fn29//77Wa7bpEkTSVJ8fHyWyzPuNT5//rxZ+fnz53O8D9nX19c003XGCwAAAMgvVwzDGVy57YClHB6I75Wenm52P+/dDhw4IEkqX758lsurVaum4OBgbd682VSWnJysnTt3mt2XDAAAABQ0VwuUYTWLZipztWMA8sqhQ6ZHjx6tDh06qHLlyrp69apWrlypmJgYbdy4UcePH9fKlSvVsWNHlSpVSnFxcRo6dKhatGih8PBwUx1hYWGaNm2annjiCRkMBr322muaMmWKatSoYXrsUoUKFdS1a1fHHSgAAAA8iqsFyYwwnPHfw0dTTMsyjoUh1HBHDg3EFy5cUHR0tM6dO6egoCCFh4dr48aNatu2rU6dOqXvv/9ec+bM0bVr1xQSEqInn3xSY8aMMavjyJEjSkpKMr0fMWKErl27pkGDBunKlStq3ry5NmzYwDOIAQAAUOBcLQhLWfcMh9UsahaKJe4rhnsyGI1Go6Mb4WySk5MVFBSkdtE75e2T+QMCAAAAuJerheGsgvC97g3FkmN6im/fStF3y5soKSnJpeb7ycgVs79Kkn+Afdt941qyhj0Z5HLnzN6c7h5iAAAAwNW4YxjObj1XO1YgJwRiAAAAIB9cLSBmFXJrVLhpelmyvqsdM5Adhz+HGAAAAHBVrhQMs+sVvjcEZ7w/dvZ/c/Aw2RbcFT3EAAAAQB6Fhld2+TCcXY/w3cstqceVzgNwLwIxAAAAkAeuFgCzC8OWIBTD3RGIAQAAAAu5UvALq1nU4jBcrcgZ08uS9QnFcBcEYgAAAMACrhT4LL1fWFKmEJzfUOxK5wkgEAMAAAC5cKWQl5f7hbMKv9mVZ1VHdr3QrnS+4NkIxAAAAEA2XKnH05oh0jlhCDU8AYEYAAAAyIIrBbr8DJHOUPbibyp78TeL1icUw10QiAEAAIB7uFKQy+8Q6XuDMKEYnoRADAAAANzFlQJcXh6plF0YzkpBhGJXOq/wHARiAAAA4P9zldBmi/uFswvDOS1nsi24GwIxAAAAPJ4r9WDaeoh0hsLHD6nw8UO5rstkW3AnBGIAAAB4NFcKaAU1RPruIHxvKM5uO0Kxe/roo49UtWpV+fn5qUmTJtq1a1eO68+ZM0e1atWSv7+/QkJCNHToUN28mfXPpDMiEAMAAMBjuUowK6gh0ln1CmeUW7I9odi9fP755xo2bJjGjx+vffv2qX79+oqMjNSFCxeyXH/lypUaNWqUxo8fr//+979atGiRPv/8c7355pt2brn1CMQAAADwSK4SyGz1SKV7ZRV6c1vOZFvubfbs2Ro4cKD69eunOnXqaP78+SpSpIgWL16c5frbt29Xs2bN1KtXL1WtWlXt2rVTz549c+1VdiYEYgAAAHgcVwlhBXm/8L2S9v9q0XrZheJ7989kW84hOTnZ7JWamprlerdu3dLevXvVpk0bU1mhQoXUpk0bxcbGZrlN06ZNtXfvXlMA/v3337Vu3Tp17NjR9gdSQLwc3QAAAADAXlwleNm7VzgjDCft/1VBDe7Pcpu0++pmqvtC6TqZ2nLiesVMbT521s+sLKxmUR0+mmJWFhpeWfFxCVm2z9UdjU+Rr799+yJTb/x9fkNCQszKx48frwkTJmRa/+LFi0pLS1O5cuXMysuVK6fDhw9nuY9evXrp4sWLat68uYxGo+7cuaMXXniBIdMAAACAsyEM5xyG736f397ie3FfseOcOnVKSUlJptfo0aNtVndMTIymTp2qjz/+WPv27dPXX3+t//znP5o8ebLN9lHQ6CEGAACA23OVoGWPWaQzZBV6712eVW/x3T3FGfukp9h5BQYGKjAwMNf1SpcurcKFC+v8+fNm5efPn1dwcHCW24wdO1bPPvusnnvuOUlSvXr1dO3aNQ0aNEhvvfWWChVy/v5X528hAAAAkA+E4byH4ZzWK4ie4nuPncm27M/Hx0cNGzbU5s2bTWXp6enavHmzIiIistzm+vXrmUJv4cKFJUlGo7HgGmtD9BADAADALblKoHL0EGlJ+uOnY6Z/V2leI9P6+ekplmTWW5xxXPQWO59hw4apT58+atSokRo3bqw5c+bo2rVr6tevnyQpOjpaFStW1LRp0yRJnTt31uzZs9WgQQM1adJE8fHxGjt2rDp37mwKxs6OQAwAAAC348ph2N69wneH4Yz3loZiyT6Tbf13T9bHCdt65pln9Oeff2rcuHFKTEzUAw88oA0bNpgm2kpISDDrER4zZowMBoPGjBmjM2fOqEyZMurcubPefvttRx1CnhmMrtKXbUfJyckKCgpSu+id8vbJ+i92AAAAcE6EYevD8L3uDcaSMgVjSZl6i6XMoVhSplAsZe4plpQpFKfeSNb8kRWVlJRk0f2wziIjV7zwzhn5+tu33a56zuyNe4gBAAAAAB6JQAwAAAC34ir3nN7bCypl3VsqZd2zmlUPrJR1b21WvbpS1j3AGeXO0jssSb//ejrLdgL5RSAGAACA23H1UJxVUMwuFGcVPrMLxVkF2nuDb3ZB+N5t0+6ra5cw7CrXEq6JQAwAAAC3FB+X4BJhKqsQKGUdGLMKllLWITSrsCpl3cub0SOc317he9tx4npFwjCcGoEYAAAAbs0VQtXhoykWD6HOLmRmF4rzMoTakvXy2yt87zFld+yucN3g+gjEAAAAcHuuEq6cYQj13cstqacghki7yvWC6yMQAwAAwCO4SshyhiHUjgzDgD0RiAEAAOAxXCVwOWoINZNnwdMQiAEAAOBRXGlIrqMfzZSfybOyG+pNGIYzIRADAADAI7lKCLPnfcW5Lc9vr/C9x+JKf5yAeyIQAwAAwGO5Shiz9xBqhkjDUxCIAQAA4NFcKZjZYwg1YRiehEAMAAAAj+dKAa0gh1DfG4azW5cwDHdBIAYAAADkWvezFtSjmXJbzuRZcDcEYgAAAOAurhLaCuq+4uzKmTwL7ohADAAAANzDlcKbrYdQM0QanoRADAAAAGTBlYKcrYZQE4bhaQjEAAAAQDZcaaivLYZQZ7WeJfURhuGqCMQAAABALlwp3OX30UwZ5UyeBU9AIAYAAAAs4EohLz/3FTN5FjwJgRgAAACwkCuFPWuGUDNEGp6GQAwAAADkgasFv7wMobZkPcIw3AmBGAAAAMgjVxsanJch1Hcvt6QeVzoPwL0IxAAAAICVXCkMWjqEOqugnN22rnT8QFa8HN0AAAAAwJXFxyUoNLyyo5thscNHUxRWs6hZWU49xQRhuDN6iAEAAIB8crWAmFXItXQ9VztWICcEYgAAAMAGXC0o5haKCcPwBARiAAAAwEZccbKtrIIvYRiegkAMAAAA2JirhceMAMzkWfA0BGIAAACgALhaiMwuCLvacQB5QSAGAAAACogrh0lXbjtgKQIxAAAAUIBcMVi6YpsBaxCIAQAAgALmSkOPXaWdgC14OboBAAAAgKcgbALOhR5iAAAAAIBHIhADAAAAADwSgRgAAAAA4JEIxAAAAAAAj0QgBgAAAAB4JAIxAAAAAMAjEYgBAAAAAB6JQAwAAAAA8EgEYgAAAACARyIQAwAAAAA8EoEYAAAAAOCRHBqI582bp/DwcAUGBiowMFARERFav369JOny5ct6+eWXVatWLfn7+6ty5cp65ZVXlJSUlGOdffv2lcFgMHu1b9/eHocDAAAAAHAhXo7ceaVKlTR9+nTVqFFDRqNRy5YtU5cuXbR//34ZjUadPXtW7777rurUqaM//vhDL7zwgs6ePavVq1fnWG/79u21ZMkS03tfX9+CPhQAAAAAgItxaCDu3Lmz2fu3335b8+bN044dOzRgwAB99dVXpmX33Xef3n77bfXu3Vt37tyRl1f2Tff19VVwcHCBtRsAAAAA4Pqc5h7itLQ0rVq1SteuXVNERESW6yQlJSkwMDDHMCxJMTExKlu2rGrVqqXBgwfr0qVLOa6fmpqq5ORksxcAAAAAwL05tIdYkg4ePKiIiAjdvHlTRYsW1Zo1a1SnTp1M6128eFGTJ0/WoEGDcqyvffv26tatm6pVq6bjx4/rzTffVIcOHRQbG6vChQtnuc20adM0ceJEmxwPAAAAAMA1GIxGo9GRDbh165YSEhKUlJSk1atX65///Ke2bt1qFoqTk5PVtm1blSxZUt988428vb0trv/333/Xfffdp++//16tW7fOcp3U1FSlpqaa7S8kJETtonfK26eo9QcHAAAAIN9u30rRd8ubmEaMuork5GQFBQXphXfOyNffvu1OvZGs+SMrutw5szeHD5n28fFRaGioGjZsqGnTpql+/fp6//33TcuvXr2q9u3bq1ixYlqzZk2ewrAkVa9eXaVLl1Z8fHy26/j6+ppmus54AQAAAADcm8MD8b3S09NNvbXJyclq166dfHx89M0338jPzy/P9Z0+fVqXLl1S+fLlbd1UAAAAAIALc2ggHj16tLZt26aTJ0/q4MGDGj16tGJiYhQVFWUKw9euXdOiRYuUnJysxMREJSYmKi0tzVRHWFiY1qxZI0lKSUnRG2+8oR07dujkyZPavHmzunTpotDQUEVGRjrqMAEAAAAATsihk2pduHBB0dHROnfunIKCghQeHq6NGzeqbdu2iomJ0c6dOyVJoaGhZtudOHFCVatWlSQdOXJESUlJkqTChQsrLi5Oy5Yt05UrV1ShQgW1a9dOkydP5lnEAAAAAAAzDg3EixYtynZZq1atZMl8X3ev4+/vr40bN9qkbQAAAAAA9+Z09xADAAAAAGAPBGIAAAAAgCTpo48+UtWqVeXn56cmTZpo165dOa7/5ZdfKiwsTH5+fqpXr57WrVtnp5baBoEYAAAAAKDPP/9cw4YN0/jx47Vv3z7Vr19fkZGRunDhQpbrb9++XT179tSAAQO0f/9+de3aVV27dtWhQ4fs3HLrEYgBAAAAAJo9e7YGDhyofv36qU6dOpo/f76KFCmixYsXZ7n++++/r/bt2+uNN95Q7dq1NXnyZD344IOaO3eunVtuPQIxAAAAALip5ORks1dqamqW6926dUt79+5VmzZtTGWFChVSmzZtFBsbm+U2sbGxZutLUmRkZLbrOyOHzjINAAAAAO7u919Py9unqF33eftWiiQpJCTErHz8+PGaMGFCpvUvXryotLQ0lStXzqy8XLlyOnz4cJb7SExMzHL9xMTEfLTcvgjEAAAAAOCmTp06pcDAQNN7X19fB7bG+RCIAQAAAMBNBQYGmgXi7JQuXVqFCxfW+fPnzcrPnz+v4ODgLLcJDg7O0/rOiHuIAQAAnFRoeGVHNwGAh/Dx8VHDhg21efNmU1l6ero2b96siIiILLeJiIgwW1+SNm3alO36zogeYgAAACdzdxAODa+s+LgEB7YGgKcYNmyY+vTpo0aNGqlx48aaM2eOrl27pn79+kmSoqOjVbFiRU2bNk2S9Oqrr6ply5aaNWuWOnXqpFWrVmnPnj1asGCBIw8jTwjEAAAATiSrXmFCMQB7eOaZZ/Tnn39q3LhxSkxM1AMPPKANGzaYJs5KSEhQoUL/G2TctGlTrVy5UmPGjNGbb76pGjVqaO3atapbt66jDiHPDEaj0ejoRjib5ORkBQUFqV30TrvPBgcAADxXbkOkCcXwVLdvpei75U2UlJRk0f2wzsKRucJVz5m9cQ8xAACAE7DkfmHuKQYA2yIQAwAAOFhegm5oeGWCMQDYCIEYAADAQfITbgnFAJB/BGIAAAAHsEWgJRQDQP4QiAEAAOzMlkGWUAwA1iMQAwAA2FFBBFjuKwYA6xCIAQAA7MAeoZVQDAB5QyAGAAAoYPYMqoRiALAcgRgAAKAAOSKgEooBwDIEYgAAgALiyGBKKAaA3BGIAQAACoAzBFIm2wKAnBGIAQAAbMgZQ6iztQcAnAWBGAAAwEacOXg6c9sAwFEIxAAAADbgCoHTFdoIAPZEIAYAAMgnVwqartRWAChoBGIAAIB8cMWA6Yz3OQOAIxCIAQAArOAOodLV2w8A+UUgBgAAyCN3CpLudCwAkFcEYgAAgDxwxwDpjscEAJYgEAMAAFjInYOjOx8bAGSHQAwAAGABTwiM7nBfNADkBYEYAAAgB54YEj3teAF4LgIxAABANjw5GHrysQPwHARiAACALBAIOQcA3B+BGAAA4B4Ewf/hXABwZwRiAACAuxAAM/PE+6gBeAYCMQAAgAh9luD8AHA3BGIAAODxCHqW41wBcCcEYgAA4NEIeHnHOQPgLgjEAAAAyDNCMQB3QCAGAABAnsXHJTi6CQCQbwRiAADg0eLjEgh3ecT5AuAuCMQAAAAi5FmCPx4AcDcEYgAAgP+PsJc9zg0Ad0QgBgAAuAvBLzPOCQB3RSAGAAC4BwHwfzgXANwZgRgAACAL3C9LGAbg/gjEAAAAOfDEUMgfAwB4CgIxAABALjwpHHrSsQIAgRgAAMACnhAUPeEYAeBuBGIAAAALuXNgdOdjA4DsEIgBAADywB3vr3W34wEASxGIAQAArOAOIdIdwz0A5AWBGAAAwEquHCZdue0AYCsEYgAAgHxwxWDpim0GgIJAIAYAAMgnVxp67CrtBAB7IBADAADYiLOHTWdvHwDYG4EYAADAhpwxdLpSDzYA2BOBGAAAwMacKXw6U1sAwNkQiAEAAAqAMwRRZ2gDADgzAjEAAEABceRQZcIwAOSOQAwAAFDA7BlOuV8YACxHIAYAALADe4RUgjAA5A2BGAAAwE4KMrAShgEg7wjEAAAAdlQQwZUwDADWIRADAADYmS3v8yUMA4D1CMQAAAAOkp8wy+RZAJB/BGIAAAAHsibUEoQBwDYcGojnzZun8PBwBQYGKjAwUBEREVq/fr1p+c2bNzVkyBCVKlVKRYsW1ZNPPqnz58/nWKfRaNS4ceNUvnx5+fv7q02bNjp27FhBHwoAAIDV8hJwCcMAHO3y5cuKiopSYGCgihcvrgEDBiglJSXHbRITE/Xss88qODhYAQEBevDBB/XVV1/ZqcXZc2ggrlSpkqZPn669e/dqz549euyxx9SlSxf9+uuvkqShQ4fq3//+t7788ktt3bpVZ8+eVbdu3XKsc8aMGfrggw80f/587dy5UwEBAYqMjNTNmzftcUgAAABWsSToEoYBOIOoqCj9+uuv2rRpk7799ltt27ZNgwYNynGb6OhoHTlyRN98840OHjyobt26qXv37tq/f7+dWp01g9FoNDq0BfcoWbKkZs6cqaeeekplypTRypUr9dRTT0mSDh8+rNq1ays2NlYPP/xwpm2NRqMqVKig119/XcOHD5ckJSUlqVy5clq6dKl69OhhURuSk5MVFBSkdtE75e1T1HYHBwAAYIHQ8MqZygjD8GS3b6Xou+VNlJSUpMDAQEc3x2KOzBUFdc7++9//qk6dOtq9e7caNWokSdqwYYM6duyo06dPq0KFClluV7RoUc2bN0/PPvusqaxUqVJ655139Nxzz9msfXnl5bA93yMtLU1ffvmlrl27poiICO3du1e3b99WmzZtTOuEhYWpcuXK2QbiEydOKDEx0WyboKAgNWnSRLGxsdkG4tTUVKWmppreJyUlSZLu3Mq52x8AAKAg/HfPb6p+fyVJ0u+/nnZwawDHy/j/cifry7PYndv2zxUZ+0xOTjYr9/X1la+vr9X1xsbGqnjx4qYwLElt2rRRoUKFtHPnTj3xxBNZbte0aVN9/vnn6tSpk4oXL64vvvhCN2/eVKtWraxuiy04PBAfPHhQERERunnzpooWLao1a9aoTp06OnDggHx8fFS8eHGz9cuVK6fExMQs68ooL1eunMXbSNK0adM0ceLETOU/rGqdx6MBAAAAUFCuXr2qoKAgRzfDYj4+PgoODtYP/3JMrihatKhCQkLMysaPH68JEyZYXWdiYqLKli1rVubl5aWSJUvmmLm++OILPfPMMypVqpS8vLxUpEgRrVmzRqGhoVa3xRYcHohr1aqlAwcOKCkpSatXr1afPn20detWu7Zh9OjRGjZsmOn9lStXVKVKFSUkJLjUL5y7SU5OVkhIiE6dOuVSQ2PcDdfBOXAdHI9r4By4Ds6B6+AcPOk6GI1GXb16NdvhuM7Kz89PJ06c0K1btxyyf6PRKIPBYFaWXe/wqFGj9M477+RY33//+1+r2zJ27FhduXJF33//vUqXLq21a9eqe/fu+vHHH1WvXj2r680vhwdiHx8f018FGjZsqN27d+v999/XM888o1u3bunKlStmvcTnz59XcHBwlnVllJ8/f17ly5c32+aBBx7Itg3ZDRsICgpy+w8XV5AxCzkci+vgHLgOjsc1cA5cB+fAdXAOnnIdXLWjys/PT35+fo5uRq5ef/119e3bN8d1qlevruDgYF24cMGs/M6dO7p8+XK2Oe348eOaO3euDh06pPvvv1+SVL9+ff3444/66KOPNH/+fJscgzUcHojvlZ6ertTUVDVs2FDe3t7avHmznnzySUnSkSNHlJCQoIiIiCy3rVatmoKDg7V582ZTAE5OTtbOnTs1ePBgex0CAAAAALiUMmXKqEyZMrmuFxERoStXrmjv3r1q2LChJOmHH35Qenq6mjRpkuU2169flyQVKmT+kKPChQsrPT09ny3PH4c+dmn06NHatm2bTp48qYMHD2r06NGKiYlRVFSUgoKCNGDAAA0bNkxbtmzR3r171a9fP0VERJhNqBUWFqY1a9ZIkgwGg1577TVNmTLFNJ13dHS0KlSooK5duzroKAEAAADAPdSuXVvt27fXwIEDtWvXLv3888966aWX1KNHD9OQ9jNnzigsLEy7du2S9HdmCw0N1fPPP69du3bp+PHjmjVrljZt2uTwnObQHuILFy4oOjpa586dU1BQkMLDw7Vx40a1bdtWkvTee++pUKFCevLJJ5WamqrIyEh9/PHHZnUcOXLENCu0JI0YMULXrl3ToEGDdOXKFTVv3lwbNmzI0zAFX19fjR8/Pl+zryH/uA7OgevgHLgOjsc1cA5cB+fAdXAOXAc4yooVK/TSSy+pdevWprz2wQcfmJbfvn1bR44cMfUMe3t7a926dRo1apQ6d+6slJQUhYaGatmyZerYsaOjDkOSEz6HGAAAAAAAe3DokGkAAAAAAByFQAwAAAAA8EgEYgAAAACARyIQAwAAAAA8kkcE4rfffltNmzZVkSJFVLx48SzXeeWVV9SwYUP5+vqanmF8t5MnT8pgMGR67dixI8d9JyQkqFOnTipSpIjKli2rN954Q3fu3LHBUbkeW1yHmJgYdenSReXLl1dAQIAeeOABrVixItd9Z3XtVq1alc8jck22uA6SFBcXp0ceeUR+fn4KCQnRjBkzct03vw9/s+Qa5Hau+vbtm+XPdcbD7rNi7eeYu7LFdYiJicnynCYmJua4b2t+f9yVLa7D119/rbZt26pMmTIKDAxURESENm7cmON++X0wZ4vrIP39O/Hggw/K19dXoaGhWrp0aa775vche/v27VPbtm1VvHhxlSpVSoMGDVJKSopp+dKlS7P8OTYYDLpw4UK29VatWjXT+tOnT7fHIQFOxyMC8a1bt/T0009r8ODBOa7Xv39/PfPMMzmu8/333+vcuXOmV8bDqLOSlpamTp066datW9q+fbuWLVumpUuXaty4cVYdh6uzxXXYvn27wsPD9dVXXykuLk79+vVTdHS0vv3221z3v2TJErNr5+hnnjmKLa5DcnKy2rVrpypVqmjv3r2aOXOmJkyYoAULFmRbH78P/5PbNbDkXL3//vtmP8+nTp1SyZIl9fTTT+e6/7x8jrkzW1yHDEeOHDE7p2XLls12v9b8/rgzW1yHbdu2qW3btlq3bp327t2rRx99VJ07d9b+/ftz3T+/D3+zxXU4ceKEOnXqpEcffVQHDhzQa6+9pueeey7HP07w+5C9s2fPqk2bNgoNDdXOnTu1YcMG/frrr+rbt69pnWeeecbs5/fcuXOKjIxUy5Ytc/wckqRJkyaZbffyyy8X8BEBTsroQZYsWWIMCgrKcZ3x48cb69evn6n8xIkTRknG/fv3W7y/devWGQsVKmRMTEw0lc2bN88YGBhoTE1Ntbged5Of65CVjh07Gvv165fjOpKMa9assayBHiI/1+Hjjz82lihRwuzneOTIkcZatWplWxe/D5lldw2sOVdr1qwxGgwG48mTJ7PdnzWfY54gP9dhy5YtRknGv/76y+L9WfP74wls+ftgNBqNderUMU6cODHb5fw+ZC0/12HEiBHG+++/32y7Z555xhgZGZnt/vh9yN4nn3xiLFu2rDEtLc1UFhcXZ5RkPHbsWJbbXLhwwejt7W1cvnx5jnVXqVLF+N5779myuYDL8ogeYlt6/PHHVbZsWTVv3lzffPNNjuvGxsaqXr16KleunKksMjJSycnJ+vXXXwu6qR4jKSlJJUuWzHW9IUOGqHTp0mrcuLEWL14sI4/gtlpsbKxatGghHx8fU1lkZKSOHDmiv/76K9tt+H2wjDXnatGiRWrTpo2qVKmSa/15+RzzZHm5Dg888IDKly+vtm3b6ueff8613rz+/ngya34f0tPTdfXqVYu+G/h9sIwl1yE2NlZt2rQx2y4yMlKxsbE51svvQ9ZSU1Pl4+OjQoX+97/r/v7+kqSffvopy22WL1+uIkWK6Kmnnsq1/unTp6tUqVJq0KCBZs6c6ZG3MAGShwyZtoWiRYtq1qxZ+vLLL/Wf//xHzZs3V9euXXP88kxMTDT74pBkep/b/WWwzBdffKHdu3erX79+Oa43adIkffHFF9q0aZOefPJJvfjii/rwww/t1Er3Y83PNr8PlsvruTp79qzWr1+v5557Lsd6rfkc82SWXIfy5ctr/vz5+uqrr/TVV18pJCRErVq10r59+/JVL/7HmvP17rvvKiUlRd27d8+2Xn4f8saS65DdOsnJybpx44bV9Xqqxx57TImJiZo5c6Zu3bqlv/76S6NGjZIknTt3LsttFi1apF69epmCc3ZeeeUVrVq1Slu2bNHzzz+vqVOnasSIETY/BsAVuGwgHjVqVLaTCGS8Dh8+bLP9lS5dWsOGDVOTJk300EMPafr06erdu7dmzpxps324Intfh7tt2bJF/fr108KFC3OcSEiSxo4dq2bNmqlBgwYaOXKkRowY4VbXzpHXAX9z5DVYtmyZihcvnut98Z7wOWbv61CrVi09//zzatiwoZo2barFixeradOmeu+992y2D1fkyN+HlStXauLEifriiy9yvIeS3we+GxzF0uty//33a9myZZo1a5aKFCmi4OBgVatWTeXKlTPrNc4QGxur//73vxowYECubRg2bJhatWql8PBwvfDCC5o1a5Y+/PBDpaamFsQhA07Ny9ENsNbrr79uNqlAVqpXr16gbWjSpIk2bdqU7fLg4GDt2rXLrOz8+fOmZe7AUddh69at6ty5s9577z1FR0fnefsmTZpo8uTJSk1Nla+vr83bZ2/2vg7BwcGmn+UMuf1su/vvgy2vQV7OldFo1OLFi/Xss8+aDTm0VG6fY67GUdfhbo0bN852OGPGtnn9/XE1jroOq1at0nPPPacvv/wy09BdS/D7kD1LrkN2P9uBgYHZ9lh6wu/DvfJyXXr16qVevXrp/PnzCgj4f+3df0xVZRzH8c+9KF6uQMC4IhnaYko4p1QTdhcLc2KZbjL/6I/KcOmsP0q3WiGzSVZYtpU50cxIamVpzaTWzAbpNaJFzXGJNnBAUFAsWxslWyyIb38477wKCsblh/f92u4fh/Oc53nO+fKcs+895z5nqhwOh1555ZUB41ZaWqqMjIyrmhguKytLfX19amtrU1pa2rC3ByayCZsQezweeTyeMe2D3+9XcnLyoOu9Xq+Ki4t15syZwLfUFRUVio2N1dy5c0ermyE1FnHw+XxasWKFtm/frvXr119VHX6/X/Hx8ddEMiyNfhy8Xq82b96s3t5eTZ48WdK5/+20tDTFx8cPus21PB5GMgbDOVYnT55Uc3PzkO4IDORK57GJZqzicKGhXBuGO34mmrGIw/vvv6+HHnpIBw8e1PLly6+qLcbD4IYSB6/Xq6NHjwZtV1FRIa/Xe9l6r/XxcLGricv5x8j3798vl8ul3NzcoPXd3d364IMP9MILL1xVn/x+v5xO5xVnpgauSWM9q9do+Omnn6y2tta2bt1q0dHRVltba7W1tXb27NlAmaamJqutrbWHH37Y5syZEyhzftbDt956y9577z1raGiwhoYGKy4uNqfTafv37w/U8dFHHwXNitjX12fz5s2zpUuXmt/vt2PHjpnH47HCwsLR2/lxZCTicPz4cXO73VZYWGidnZ2Bzx9//BGo4+I4fPLJJ/bGG29YfX29NTU12Z49e8ztdtuWLVtGb+fHkZGIQ1dXlyUlJdnq1avthx9+sIMHD5rb7bbXX389UAfjYXBXisFwjtUDDzxgWVlZA7aza9cuW7x4cWB5KOexcDIScdixY4eVl5dbU1OT1dfX28aNG83pdFplZWWgzMVxGMr4CScjEYcDBw7YpEmTbPfu3UHXhq6urkAZxsPljUQcfvzxR3O73fbkk09aQ0OD7d692yIiIuzYsWOBMoyH4dm1a5edOnXKTp8+bSUlJRYVFWU7d+68pFxpaam5XK4BZ7yvqamxtLQ06+joMDOzr7/+2nbs2GF+v99aWlrs3XffNY/HYw8++GCodwcYl8IiIc7PzzdJl3xOnDgRKJOTkzNgmdbWVjM7d+FMT083t9ttsbGxlpmZaR9++GFQO2VlZXbxdwxtbW22bNkyi4qKssTERHviiSest7c31Ls8Lo1EHAarIycnJ1DHxXH47LPPLCMjw6Kjo23q1Km2YMEC27t3b9BrDMLJSMTBzKyurs6ys7NtypQpNmPGDHvxxReD2mE8DG4oMRjKserq6rKoqCjbt2/fgO0UFRXZrFmzAstDOY+Fk5GIw/bt2y01NdVcLpclJCTYokWL7Pjx40HtXBwHsyuPn3AyEnEY7JyVn58fKMN4uLyROi+dOHHCMjIyLDIy0m666SYrKysLWs94GJ7Vq1dbQkKCRUZG2vz58wd9nZLX67X77rtvwHXnXw93/hp+6tQpy8rKsuuuu85cLpelp6fbtm3brKenJ1S7AYxrDjPePQMAAAAACD8TdpZpAAAAAAD+DxJiAAAAAEBYIiEGAAAAAIQlEmIAAAAAQFgiIQYAAAAAhCUSYgAAAABAWCIhBgAAAACEJRJiAAAAAEBYIiEGAFyz2tra5HA45Pf7Q1K/w+FQeXl5SOoGAAChR0IMAAiZNWvWKC8vb8zaT0lJUWdnp+bNmydJ8vl8cjgc6urqGrM+AQCA8WPSWHcAAIBQiYiI0PTp08e6GwAAYJziDjEAYEycPHlSmZmZmjJlipKTk7Vp0yb19fUF1i9atEgbNmzQU089pYSEBE2fPl3PPPNMUB2NjY3Kzs6Wy+XS3LlzVVlZGfQY84WPTLe1tenOO++UJMXHx8vhcGjNmjWSpBtvvFGvvvpqUN0ZGRlB7TU1NemOO+4ItFVRUXHJPrW3t+vee+9VXFycEhIStHLlSrW1tf3fQwUAAEKEhBgAMOp++eUX3XPPPVq4cKHq6ur02muv6c0339Tzzz8fVO7tt9/W1KlTVVNTo5deeknPPvtsIBH9999/lZeXJ7fbrZqaGu3bt0+bN28etM2UlBQdPnxYknT69Gl1dnZq586dQ+pvf3+/Vq1apcjISNXU1Gjv3r0qKCgIKtPb26u77rpLMTExqqqqUnV1taKjo3X33Xfrn3/+Gc7hAQAAo4RHpgEAo27Pnj1KSUlRSUmJHA6Hbr75Zv36668qKCjQli1b5HSe+752/vz5KioqkiTNnj1bJSUl+uKLL5Sbm6uKigq1tLTI5/MFHosuLi5Wbm7ugG1GREQoISFBkjRt2jTFxcUNub+VlZVqbGzU559/ruuvv16StG3bNi1btixQ5tChQ+rv71dpaakcDockqaysTHFxcfL5fFq6dOnwDhIAAAg5EmIAwKhraGiQ1+sNJI6SdPvtt6u7u1sdHR2aOXOmpHMJ8YWSk5N15swZSefu8qakpAT9RjgzMzNk/U1JSQkkw5Lk9XqDytTV1am5uVkxMTFBf+/p6VFLS0tI+gUAAP4fEmIAwLg1efLkoGWHw6H+/v4Rb8fpdMrMgv7W29s7rDq6u7t122236cCBA5es83g8/6t/AAAgNEiIAQCjLj09XYcPH5aZBe4SV1dXKyYmRjfccMOQ6khLS1N7e7t+++03JSUlSZK+++67y24TGRkp6dzvjy/k8XjU2dkZWP7rr7/U2toa1N/29nZ1dnYqOTlZkvTNN98E1XHrrbfq0KFDmjZtmmJjY4e0DwAAYGwxqRYAIKT+/PNP+f3+oM/69evV3t6uxx57TI2Njfr4449VVFSkxx9/PPD74SvJzc1Vamqq8vPz9f3336u6ulpPP/20JAU9in2hWbNmyeFw6NNPP9Xvv/+u7u5uSdLixYv1zjvvqKqqSvX19crPz1dERERguyVLlmjOnDnKz89XXV2dqqqqLpnA6/7771diYqJWrlypqqoqtba2yufzacOGDero6LiaQwcAAEKMhBgAEFI+n0+33HJL0Oe5557T0aNH9e2332rBggV65JFHtHbt2kBCOxQREREqLy9Xd3e3Fi5cqHXr1gWSVJfLNeA2M2bM0NatW7Vp0yYlJSXp0UcflSQVFhYqJydHK1as0PLly5WXl6fU1NTAdk6nU0eOHNHff/+tzMxMrVu3TsXFxUF1u91uffnll5o5c6ZWrVql9PR0rV27Vj09PdwxBgBgnHLYxT+aAgBggqqurlZ2draam5uDEloAAICBkBADACasI0eOKDo6WrNnz1Zzc7M2btyo+Ph4ffXVV2PdNQAAMAEwqRYAYMI6e/asCgoK9PPPPysxMVFLlizRyy+/PNbdAgAAEwR3iAEAAAAAYYlJtQAAAAAAYYmEGAAAAAAQlkiIAQAAAABhiYQYAAAAABCWSIgBAAAAAGGJhBgAAAAAEJZIiAEAAAAAYYmEGAAAAAAQlv4D39gZ1M2D5wsAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot results\n", + "plt.figure(figsize=(12, 8))\n", + "ef_per_m = result.data[\"ef_per_m\"].isel(time=0, level=0)\n", + "plt.contourf(result.data.longitude, result.data.latitude, ef_per_m, cmap='coolwarm')\n", + "plt.colorbar(label=\"Energy Forcing per meter (W/m)\")\n", + "plt.title(\"CocipGrid Energy Forcing\")\n", + "plt.xlabel(\"Longitude\")\n", + "plt.ylabel(\"Latitude\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "MetDataset with data:

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 87kB\n",
+       "Dimensions:       (longitude: 20, latitude: 20, level: 2, time: 9)\n",
+       "Coordinates:\n",
+       "  * longitude     (longitude) float64 160B -115.0 -114.0 -113.0 ... -97.0 -96.0\n",
+       "  * latitude      (latitude) float64 160B 30.0 31.0 32.0 33.0 ... 47.0 48.0 49.0\n",
+       "  * level         (level) float64 16B 250.0 300.0\n",
+       "  * time          (time) datetime64[ns] 72B 2022-03-01 ... 2022-03-02\n",
+       "    air_pressure  (level) float32 8B 2.5e+04 3e+04\n",
+       "    altitude      (level) float32 8B 1.036e+04 9.164e+03\n",
+       "Data variables:\n",
+       "    contrail_age  (longitude, latitude, level, time) float32 29kB 0.0 ... 0.0\n",
+       "    ef_per_m      (longitude, latitude, level, time) float64 58kB 0.0 ... 0.0\n",
+       "Attributes: (12/13)\n",
+       "    humidity_scaling_name:     histogram_matching\n",
+       "    humidity_scaling_formula:  era5_quantiles -> iagos_quantiles\n",
+       "    azimuth:                   0.0\n",
+       "    segment_length:            1000.0\n",
+       "    met_source_provider:       ECMWF\n",
+       "    met_source_dataset:        ERA5\n",
+       "    ...                        ...\n",
+       "    aircraft_type:             B737\n",
+       "    description:               Gridded Contrail Cirrus Prediction Model\n",
+       "    max_age:                   20 hours\n",
+       "    dt_integration:            30 minutes\n",
+       "    pycontrails_version:       0.54.0\n",
+       "    ap_model:                  PSGrid
" + ], + "text/plain": [ + "MetDataset with data:\n", + "\n", + " Size: 87kB\n", + "Dimensions: (longitude: 20, latitude: 20, level: 2, time: 9)\n", + "Coordinates:\n", + " * longitude (longitude) float64 160B -115.0 -114.0 -113.0 ... -97.0 -96.0\n", + " * latitude (latitude) float64 160B 30.0 31.0 32.0 33.0 ... 47.0 48.0 49.0\n", + " * level (level) float64 16B 250.0 300.0\n", + " * time (time) datetime64[ns] 72B 2022-03-01 ... 2022-03-02\n", + " air_pressure (level) float32 8B 2.5e+04 3e+04\n", + " altitude (level) float32 8B 1.036e+04 9.164e+03\n", + "Data variables:\n", + " contrail_age (longitude, latitude, level, time) float32 29kB 0.0 ... 0.0\n", + " ef_per_m (longitude, latitude, level, time) float64 58kB 0.0 ... 0.0\n", + "Attributes: (12/13)\n", + " humidity_scaling_name: histogram_matching\n", + " humidity_scaling_formula: era5_quantiles -> iagos_quantiles\n", + " azimuth: 0.0\n", + " segment_length: 1000.0\n", + " met_source_provider: ECMWF\n", + " met_source_dataset: ERA5\n", + " ... ...\n", + " aircraft_type: B737\n", + " description: Gridded Contrail Cirrus Prediction Model\n", + " max_age: 20 hours\n", + " dt_integration: 30 minutes\n", + " pycontrails_version: 0.54.0\n", + " ap_model: PSGrid" + ] + }, + "execution_count": 42, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean Energy Forcing per meter: 3626066.1523160664\n", + "Max Energy Forcing per meter: 615089695.7354624\n", + "Min Energy Forcing per meter: -16682450.716573667\n", + "Mean Contrail Age: 0.1274999976158142\n" + ] + } + ], + "source": [ + "# work with the results\n", + "print(\"Mean Energy Forcing per meter:\", ef_per_m.mean().item())\n", + "print(\"Max Energy Forcing per meter:\", ef_per_m.max().item())\n", + "print(\"Min Energy Forcing per meter:\", ef_per_m.min().item())\n", + "\n", + "contrail_age = result.data[\"contrail_age\"].isel(time=0, level=0)\n", + "print(\"Mean Contrail Age:\", contrail_age.mean().item())" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "venv", + "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.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 1d2eb62247fc7d2b48dbd421c517e534fe3716ac Mon Sep 17 00:00:00 2001 From: Ellie Daw Date: Tue, 1 Oct 2024 12:04:14 -0500 Subject: [PATCH 2/3] Clears warning output from cell --- docs/notebooks/CoCiPGrid.ipynb | 24 +++++------------------- 1 file changed, 5 insertions(+), 19 deletions(-) diff --git a/docs/notebooks/CoCiPGrid.ipynb b/docs/notebooks/CoCiPGrid.ipynb index d92de6243..2d0ee0f35 100644 --- a/docs/notebooks/CoCiPGrid.ipynb +++ b/docs/notebooks/CoCiPGrid.ipynb @@ -73,23 +73,9 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/ellie/devspace/breakthrough_analysis/venv/lib/python3.10/site-packages/pycontrails/models/cocipgrid/cocip_grid.py:2426: UserWarning: Met data 'rad' does not cover the source domain along the time axis. This causes some interpolated values to be nan, leading to meaningless results.\n", - " warnings.warn(\n", - "/Users/ellie/devspace/breakthrough_analysis/venv/lib/python3.10/site-packages/pycontrails/models/cocipgrid/cocip_grid.py:2568: UserWarning: End time of parameter 'met' (2022-03-02 00:00:00) is before model end time (2022-03-02 20:30:00). Include additional time at the end of 'met' or reduce 'max_age' parameter.\n", - " warnings.warn(\n", - "/Users/ellie/devspace/breakthrough_analysis/venv/lib/python3.10/site-packages/pycontrails/models/cocipgrid/cocip_grid.py:2568: UserWarning: End time of parameter 'rad' (2022-03-01 23:30:00) is before model end time (2022-03-02 20:30:00). Include additional time at the end of 'rad' or reduce 'max_age' parameter. Note: differencing reduces time coverage when providing accumulated radiative fluxes.\n", - " warnings.warn(\n", - "CocipGrid eval: 73%|███████▎ | 65/89 [00:04<00:01, 15.73it/s]\n" - ] - } - ], + "outputs": [], "source": [ "# Create a grid source\n", "grid_source = MetDataset.from_coords(**grid_params)\n", @@ -695,7 +681,7 @@ ], "metadata": { "kernelspec": { - "display_name": "venv", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -709,9 +695,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.13" + "version": "3.11.6" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } From 032dfd8414c32e92db1ccd622640ae695438fca8 Mon Sep 17 00:00:00 2001 From: Ellie Daw Date: Tue, 1 Oct 2024 13:27:15 -0500 Subject: [PATCH 3/3] Linter :) --- docs/notebooks/CoCiPGrid.ipynb | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/notebooks/CoCiPGrid.ipynb b/docs/notebooks/CoCiPGrid.ipynb index 2d0ee0f35..0a94ab6b5 100644 --- a/docs/notebooks/CoCiPGrid.ipynb +++ b/docs/notebooks/CoCiPGrid.ipynb @@ -6,14 +6,15 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "from pycontrails.models.cocipgrid import CocipGrid\n", - "from pycontrails.datalib.ecmwf import ERA5\n", - "from pycontrails.core import MetDataset\n", + "import numpy as np\n", "import pandas as pd\n", - "from pycontrails.models.ps_model import PSGrid\n", - "from pycontrails.models.humidity_scaling import HistogramMatching" + "\n", + "from pycontrails.core import MetDataset\n", + "from pycontrails.datalib.ecmwf import ERA5\n", + "from pycontrails.models.cocipgrid import CocipGrid\n", + "from pycontrails.models.humidity_scaling import HistogramMatching\n", + "from pycontrails.models.ps_model import PSGrid" ] }, { @@ -104,7 +105,7 @@ "# Plot results\n", "plt.figure(figsize=(12, 8))\n", "ef_per_m = result.data[\"ef_per_m\"].isel(time=0, level=0)\n", - "plt.contourf(result.data.longitude, result.data.latitude, ef_per_m, cmap='coolwarm')\n", + "plt.contourf(result.data.longitude, result.data.latitude, ef_per_m, cmap=\"coolwarm\")\n", "plt.colorbar(label=\"Energy Forcing per meter (W/m)\")\n", "plt.title(\"CocipGrid Energy Forcing\")\n", "plt.xlabel(\"Longitude\")\n", @@ -694,8 +695,7 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.11.6" + "pygments_lexer": "ipython3" } }, "nbformat": 4,