{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Grid indexing\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Parcels supports `Fields` on any curvilinear `Grid`. For velocity `Fields` (`U`, `V` and `W`), there are some additional restrictions if the `Grid` is discretized as an Arakawa B- or Arakawa C-grid. That is because under the hood, Parcels uses a specific interpolation scheme for each of these grid types. This is described in [Section 2 of Delandmeter and Van Sebille (2019)](https://www.geosci-model-dev.net/12/3571/2019/gmd-12-3571-2019.html) and summarized below.\n", "\n", "The summary is: \n", "\n", "> For Arakawa B-and C-grids, Parcels requires the locations of the _corner_ points (f-points) of the grid cells for the `dimensions` dictionary of velocity `Fields`.\n", " \n", "In other words, on an Arakawa C-grid, the `[k, j, i]` node of `U` will _not_ be located at position `[k, j, i]` of `U.grid`.\n", "\n", "## Barycentric coordinates and indexing in Parcels\n", "\n", "### arakawa C-grids\n", "\n", "In a regular grid (also called an Arakawa A-grid), the velocities (`U`, `V` and `W`) and tracers (e.g. temperature) are all located in the center of each cell. But hydrodynamic model data is often supplied on staggered grids (e.g. for NEMO, POP and MITgcm), where `U`, `V` and `W` are shifted with respect to the cell centers.\n", "\n", "To keep it simple, let's take the case of a 2D Arakawa C-grid. Zonal (`U`) velocities are located at the east and west faces of each cell and meridional (`V`) velocities at the north and south. The following diagram shows a comparison of 4x4 A- and C-grids.\n", "\n", "![Arakawa Grid layouts](images/grid_comparison.png)\n", "\n", "Note that different models use different indices to relate the location of the staggered velocities to the cell centers. The default C-grid indexing notation in Parcels is the same as for **NEMO** (middle panel): `U[j, i]` is located between the cell corners `[j-1, i]` and `[j, i]`, and `V[j, i]` is located between cell corners `[j, i-1]` and `[j, i]`.\n", "\n", "Now, as you've noticed on the grid illustrated on the figure, there are 4x4 cells. The grid providing the cell corners is a 5x5 grid, but there are only 4x5 U nodes and 5x4 V nodes, since the grids are staggered. This implies that the first row of `U` data and the first column of `V` data is never used (and do not physically exist), but the `U` and `V` fields are correctly provided on a 5x5 table. If your original data are provided on a 4x5 `U` grid and a 5x4 `V` grid, you need to regrid your table to follow Parcels notation before creating a FieldSet!\n", "\n", "**MITgcm** uses a different indexing: `U[j, i]` is located between the cell corners `[j, i]` and `[j+1, i]`, and `V[j, i]` is located between cell corners `[j, i]` and `[j, i+1]`. If you use [xmitgcm](https://xmitgcm.readthedocs.io/en/latest/) to write your data to a NetCDF file, `U` and `V` will have the same dimensions as your grid (in the above case 4x4 rather than 5x5 as in NEMO), meaning that the last column of `U` and the last row of `V` are omitted. In MITgcm, these either correspond to a solid boundary, or in the case of a periodic domain, to the first column and row respectively. In the latter case, and assuming your model domain is periodic, you can use the `FieldSet.add_periodic_halo` method to make sure particles can be correctly interpolated in the last zonal/meridional cells.\n", "\n", "Parcels can take care of loading C-grid data for you, through the general `FieldSet.from_c_grid_dataset` method (which takes a `gridindexingtype` argument), and the model-specific methods `FieldSet.from_nemo` and `FieldSet.from_mitgcm`.\n", "\n", "The only information that Parcels needs is the location of the 4 corners of the cell, which are called the $f$-node. Those are the ones provided by `U.grid` (and are equal to the ones in `V.grid`). Parcels does not need the exact location of the `U` and `V` nodes, because in C-grids, `U` and `V` nodes are always located in the same position relative to the $f$-node.\n", "\n", "Importantly, the interpolation of velocities on Arakawa C-grids is done in barycentric coordinates: $\\xi$, $\\eta$ and $\\zeta$ instead of longitude, latitude and depth. These coordinates are always between 0 and 1, measured from the corner of the grid cell. This is more accurate than simple linear interpolation on curvilinear grids.\n", "\n", "When calling `FieldSet.from_c_grid_dataset()` or a method that is based on it like `FieldSet.from_nemo()`, Parcels automatically sets the interpolation method for the `U`, `V` and `W` `Fields` to `cgrid_velocity`. With this interpolation method, the velocity interpolation follows the description in [Section 2.1.2 of Delandmeter and Van Sebille (2019)](https://www.geosci-model-dev.net/12/3571/2019/gmd-12-3571-2019.html).\n", "\n", "#### NEMO Example\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from datetime import timedelta as delta\n", "from glob import glob\n", "from os import path\n", "\n", "import numpy as np\n", "\n", "from parcels import FieldSet, download_example_dataset" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Let's see how this works. We'll use the NemoNorthSeaORCA025-N006_data data, which is on an arakawa C-grid. If we create the `FieldSet` using the coordinates in the U, V and W files, we get an Error as seen below:\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [ "raises-exception" ] }, "outputs": [ { "ename": "ValueError", "evalue": "On a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U and V. See also ../examples/documentation_indexing.ipynb", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m/Users/erik/Codes/ParcelsCode/docs/examples/documentation_indexing.ipynb Cell 5\u001b[0m line \u001b[0;36m2\n\u001b[1;32m 7\u001b[0m variables \u001b[39m=\u001b[39m {\u001b[39m\"\u001b[39m\u001b[39mU\u001b[39m\u001b[39m\"\u001b[39m: \u001b[39m\"\u001b[39m\u001b[39muo\u001b[39m\u001b[39m\"\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39mV\u001b[39m\u001b[39m\"\u001b[39m: \u001b[39m\"\u001b[39m\u001b[39mvo\u001b[39m\u001b[39m\"\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39mW\u001b[39m\u001b[39m\"\u001b[39m: \u001b[39m\"\u001b[39m\u001b[39mwo\u001b[39m\u001b[39m\"\u001b[39m}\n\u001b[1;32m 8\u001b[0m dimensions \u001b[39m=\u001b[39m {\n\u001b[1;32m 9\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mU\u001b[39m\u001b[39m\"\u001b[39m: {\n\u001b[1;32m 10\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mlon\u001b[39m\u001b[39m\"\u001b[39m: \u001b[39m\"\u001b[39m\u001b[39mnav_lon\u001b[39m\u001b[39m\"\u001b[39m,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 26\u001b[0m },\n\u001b[1;32m 27\u001b[0m }\n\u001b[0;32m---> 29\u001b[0m fieldset \u001b[39m=\u001b[39m FieldSet\u001b[39m.\u001b[39;49mfrom_nemo(filenames, variables, dimensions)\n", "File \u001b[0;32m~/Codes/ParcelsCode/parcels/fieldset.py:545\u001b[0m, in \u001b[0;36mFieldSet.from_nemo\u001b[0;34m(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, chunksize, **kwargs)\u001b[0m\n\u001b[1;32m 543\u001b[0m \u001b[39mif\u001b[39;00m kwargs\u001b[39m.\u001b[39mpop(\u001b[39m'\u001b[39m\u001b[39mgridindexingtype\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mnemo\u001b[39m\u001b[39m'\u001b[39m) \u001b[39m!=\u001b[39m \u001b[39m'\u001b[39m\u001b[39mnemo\u001b[39m\u001b[39m'\u001b[39m:\n\u001b[1;32m 544\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mgridindexingtype must be \u001b[39m\u001b[39m'\u001b[39m\u001b[39mnemo\u001b[39m\u001b[39m'\u001b[39m\u001b[39m in FieldSet.from_nemo(). Use FieldSet.from_c_grid_dataset otherwise\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m--> 545\u001b[0m fieldset \u001b[39m=\u001b[39m \u001b[39mcls\u001b[39;49m\u001b[39m.\u001b[39;49mfrom_c_grid_dataset(filenames, variables, dimensions, mesh\u001b[39m=\u001b[39;49mmesh, indices\u001b[39m=\u001b[39;49mindices, time_periodic\u001b[39m=\u001b[39;49mtime_periodic,\n\u001b[1;32m 546\u001b[0m allow_time_extrapolation\u001b[39m=\u001b[39;49mallow_time_extrapolation, tracer_interp_method\u001b[39m=\u001b[39;49mtracer_interp_method,\n\u001b[1;32m 547\u001b[0m chunksize\u001b[39m=\u001b[39;49mchunksize, gridindexingtype\u001b[39m=\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mnemo\u001b[39;49m\u001b[39m'\u001b[39;49m, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m 548\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mhasattr\u001b[39m(fieldset, \u001b[39m'\u001b[39m\u001b[39mW\u001b[39m\u001b[39m'\u001b[39m):\n\u001b[1;32m 549\u001b[0m fieldset\u001b[39m.\u001b[39mW\u001b[39m.\u001b[39mset_scaling_factor(\u001b[39m-\u001b[39m\u001b[39m1.\u001b[39m)\n", "File \u001b[0;32m~/Codes/ParcelsCode/parcels/fieldset.py:656\u001b[0m, in \u001b[0;36mFieldSet.from_c_grid_dataset\u001b[0;34m(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, gridindexingtype, chunksize, **kwargs)\u001b[0m\n\u001b[1;32m 584\u001b[0m \u001b[39m\u001b[39m\u001b[39m\"\"\"Initialises FieldSet object from NetCDF files of Curvilinear NEMO fields.\u001b[39;00m\n\u001b[1;32m 585\u001b[0m \n\u001b[1;32m 586\u001b[0m \u001b[39mSee `here <../examples/documentation_indexing.ipynb>`__\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 653\u001b[0m \u001b[39m Keyword arguments passed to the :func:`Fieldset.from_netcdf` constructor.\u001b[39;00m\n\u001b[1;32m 654\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 655\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39m'\u001b[39m\u001b[39mU\u001b[39m\u001b[39m'\u001b[39m \u001b[39min\u001b[39;00m dimensions \u001b[39mand\u001b[39;00m \u001b[39m'\u001b[39m\u001b[39mV\u001b[39m\u001b[39m'\u001b[39m \u001b[39min\u001b[39;00m dimensions \u001b[39mand\u001b[39;00m dimensions[\u001b[39m'\u001b[39m\u001b[39mU\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m!=\u001b[39m dimensions[\u001b[39m'\u001b[39m\u001b[39mV\u001b[39m\u001b[39m'\u001b[39m]:\n\u001b[0;32m--> 656\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mOn a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U and V. \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 657\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mSee also ../examples/documentation_indexing.ipynb\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m 658\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39m'\u001b[39m\u001b[39mU\u001b[39m\u001b[39m'\u001b[39m \u001b[39min\u001b[39;00m dimensions \u001b[39mand\u001b[39;00m \u001b[39m'\u001b[39m\u001b[39mW\u001b[39m\u001b[39m'\u001b[39m \u001b[39min\u001b[39;00m dimensions \u001b[39mand\u001b[39;00m dimensions[\u001b[39m'\u001b[39m\u001b[39mU\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m!=\u001b[39m dimensions[\u001b[39m'\u001b[39m\u001b[39mW\u001b[39m\u001b[39m'\u001b[39m]:\n\u001b[1;32m 659\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mOn a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U, V and W. \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 660\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mSee also ../examples/documentation_indexing.ipynb\u001b[39m\u001b[39m\"\u001b[39m)\n", "\u001b[0;31mValueError\u001b[0m: On a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U and V. See also ../examples/documentation_indexing.ipynb" ] } ], "source": [ "example_dataset_folder = download_example_dataset(\"NemoNorthSeaORCA025-N006_data\")\n", "ufiles = sorted(glob(f\"{example_dataset_folder}/ORCA*U.nc\"))\n", "vfiles = sorted(glob(f\"{example_dataset_folder}/ORCA*V.nc\"))\n", "wfiles = sorted(glob(f\"{example_dataset_folder}/ORCA*W.nc\"))\n", "\n", "filenames = {\"U\": ufiles, \"V\": vfiles, \"W\": wfiles}\n", "variables = {\"U\": \"uo\", \"V\": \"vo\", \"W\": \"wo\"}\n", "dimensions = {\n", " \"U\": {\n", " \"lon\": \"nav_lon\",\n", " \"lat\": \"nav_lat\",\n", " \"depth\": \"depthu\",\n", " \"time\": \"time_counter\",\n", " },\n", " \"V\": {\n", " \"lon\": \"nav_lon\",\n", " \"lat\": \"nav_lat\",\n", " \"depth\": \"depthv\",\n", " \"time\": \"time_counter\",\n", " },\n", " \"W\": {\n", " \"lon\": \"nav_lon\",\n", " \"lat\": \"nav_lat\",\n", " \"depth\": \"depthw\",\n", " \"time\": \"time_counter\",\n", " },\n", "}\n", "\n", "fieldset = FieldSet.from_nemo(filenames, variables, dimensions)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can still load the data this way, if we use the `FieldSet.from_netcdf()` method. But because this method assumes the data is on an Arakawa-A grid, **this will mean wrong interpolation**.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "fieldsetA = FieldSet.from_netcdf(filenames, variables, dimensions)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Instead, we need to provide the coordinates of the $f$-points. In NEMO, these are called `glamf`, `gphif` and `depthw` (in MITgcm, these would be called `XG`, `YG` and `Zl`). The first two are in the `coordinates.nc` file, the last one is in the `W` file.\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "mesh_mask = f\"{example_dataset_folder}/coordinates.nc\"\n", "\n", "filenames = {\n", " \"U\": {\"lon\": mesh_mask, \"lat\": mesh_mask, \"depth\": wfiles[0], \"data\": ufiles},\n", " \"V\": {\"lon\": mesh_mask, \"lat\": mesh_mask, \"depth\": wfiles[0], \"data\": vfiles},\n", " \"W\": {\"lon\": mesh_mask, \"lat\": mesh_mask, \"depth\": wfiles[0], \"data\": wfiles},\n", "}\n", "\n", "# Note that all variables need the same dimensions in a C-Grid\n", "c_grid_dimensions = {\n", " \"lon\": \"glamf\",\n", " \"lat\": \"gphif\",\n", " \"depth\": \"depthw\",\n", " \"time\": \"time_counter\",\n", "}\n", "dimensions = {\n", " \"U\": c_grid_dimensions,\n", " \"V\": c_grid_dimensions,\n", " \"W\": c_grid_dimensions,\n", "}\n", "\n", "fieldsetC = FieldSet.from_nemo(\n", " filenames, variables, dimensions, netcdf_decodewarning=False\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Note by the way, that we used `netcdf_decodewarning=False` in the `FieldSet.from_nemo()` call above. This is to silence an expected warning because the time dimension in the `coordinates.nc` file can't be decoded by `xarray`.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the different grid points to see that indeed, the longitude and latitude of `fieldsetA.U` and `fieldsetA.V` are different.\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxMAAAGdCAYAAABpfy7RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABINklEQVR4nO3de1xUdf7H8fcwKAgC5oWYhAYVFRWVVirN3VLS1DU3zTUr21YjxWJLLcvMzFtFuVZmV/E3anbD3NTdLSvb0jR1k1DMLS9goFCUmxfwigrn9wfLJHI/zjBAr+fjMY9pzvVzvg3Oec/3e+ZYDMMwBAAAAAA15OXpAgAAAADUT4QJAAAAAKYQJgAAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJhCmAAAAABgirenC3CVoqIi/fDDDwoICJDFYvF0OQAAoBoMw9CxY8d02WWXycuL7ziB+qbBhIkffvhBYWFhni4DAACYkJ2drdDQUE+XAaCGGkyYCAgIkFT8j1FgYKCHqwEAANWRn5+vsLAw5+c4gPqlwYSJkqFNgYGBhAkAAOoZhigD9RODEwEAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJhCmAAAAABgCmECAAAAgCmECQAAAACmECYAAAAAmEKYAAAAAGAKYaIGcnKkdeuKnwEAAIBfO8JENTkckt0uxcYWPzscnq6o4crJydG6deuUQ2pzG4IxAABwBcJENeTkSOPGSUVFxa+LiqT4eE7E3MHhcMhutys2NlZ2u10OUpvLEYxrD8EYANDQESaqIT39lyBRorBQysjwTD0NVU5OjsaNG6ei/zV2UVGR4uPjORFzIYJx7SEY1x562tyPYAygIoSJamjfXvK6oKWsVikiwjP1NFTp6enOIFGisLBQGaQ2lyEY1w6Cce2hp839CMYAKkOYqIbQUCkpqThASMXPCxcWT4frtG/fXl4XpDar1aoIUpvLEIxrB8G4dtDT5n4EYwBVIUxUU1yclJVV3JWelVX8Gq4VGhqqpKQkWf+X2qxWqxYuXKhQUpvLEIxrB8G4dtDT5n4EYwBVsRiGYXi6CFfIz89XUFCQ8vLyFBgY6OlycBFycnKUkZGhiIgIgoSb5OQUn3BFRBAk3MXhcCg+Pl6FhYXOYBzHtxAulZNTPLTp/HNdq7X4Cx/e166Rk5Mju91eKlBYrVZlZWW57N9nPr+B+o0wAQBuQjB2P4ejeGhTYeEvPW1kNtdydzDm8xuo3wgTAIB6jZ4293NnMObzG6jfvD1dAAAAFyM0lBDhbqGhofSuASgXF2ADAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTLipMJCYmymKxaOLEic5po0ePlsViKfXo2bNnpdvp06dPmXUsFosGDx58MeUBAAAAcCNvsyumpKQoKSlJ3bp1KzNv4MCBWrJkifN148aNK93WypUrdebMGefrQ4cOqXv37hoxYoTZ8gAAAAC4makwcfz4cY0aNUqLFi3SE088UWa+j4+PQkJCqr295s2bl3qdnJwsPz8/wgQAAABQh5ka5pSQkKDBgwerX79+5c5fv369goOD1aFDB40dO1YHDx6s0fYdDoduvfVW+fv7V7hMQUGB8vPzSz0AAAAA1J4a90wkJydr27ZtSklJKXf+oEGDNGLECNntdmVmZmr69OmKjY1VamqqfHx8qtz+1q1b9Z///EcOh6PS5RITEzVr1qyalg8AAADARSyGYRjVXTg7O1sxMTFau3atunfvLqn44uno6GjNnz+/3HVyc3Nlt9uVnJysm2++ucp9xMfHa/Pmzdq5c2elyxUUFKigoMD5Oj8/X2FhYcrLy1NgYGB1DwkAAHhQfn6+goKC+PwG6qka9Uykpqbq4MGD6tGjh3NaYWGhNmzYoJdeekkFBQWyWq2l1rHZbLLb7UpPT69y+ydPnlRycrJmz55d5bI+Pj7V6ukAAAAA4B41ChPXX399mR6DMWPGKDIyUlOmTCkTJKTiX2bKzs6WzWarcvvvvvuuCgoKdMcdd9SkLAAAAAAeUKMwERAQoKioqFLT/P391aJFC0VFRen48eOaOXOmhg8fLpvNpqysLD366KNq2bKlhg0b5lznzjvvVOvWrZWYmFhqWw6HQ0OHDlWLFi0u4pAAAAAA1AbT95koj9Vq1c6dO7Vs2TIdPXpUNptNffv21fLlyxUQEOBc7sCBA/LyKv1DUnv37tUXX3yhtWvXurIkAAAAAG5Sowuw6zIu4AIAoP7h8xuo30zdZwIAAAAACBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMoG7JyZHWrSt+BgAAQJ1GmEDd4XBIdrsUG1v87HB4uiIAAABUgjCBuiEnRxo3TioqKn5dVCTFx9NDgfqLXjYAwK8AYQJ1Q3r6L0GiRGGhlJHhmXqAi0EvGxoSgjGAShAmUDe0by95XfB2tFqliAjP1AOYRS8bGhKCMYAqECZQN4SGSklJxQFCKn5euLB4OlCf0MuGhoJgDKAavD1dAOAUFycNGFB80hURQZBA/VTSy3Z+oKCXDfVRZcGYf58B/A89E6hbQkOlPn34oEL9RS8bGgqGnwKoBsIEALhaXJyUlVV80WpWVvFroL4hGAOoBothGIani3CF/Px8BQUFKS8vT4GBgZ4uBwCAhiEnx63DT/n8Buo3rpkAAAAVCw2lNwJAhRjmBAAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMCUiwoTiYmJslgsmjhxonPa6NGjZbFYSj169uxZ5baOHj2qhIQE2Ww2+fr6qlOnTlqzZs3FlAcAAADAjbzNrpiSkqKkpCR169atzLyBAwdqyZIlzteNGzeudFtnzpxR//79FRwcrL/97W8KDQ1Vdna2AgICzJYHAAAAwM1MhYnjx49r1KhRWrRokZ544oky8318fBQSElLt7S1evFiHDx/W5s2b1ahRI0mS3W43UxoAAACAWmJqmFNCQoIGDx6sfv36lTt//fr1Cg4OVocOHTR27FgdPHiw0u394x//UK9evZSQkKBLL71UUVFReuqpp1RYWFjhOgUFBcrPzy/1AAAAAFB7atwzkZycrG3btiklJaXc+YMGDdKIESNkt9uVmZmp6dOnKzY2VqmpqfLx8Sl3ne+++06fffaZRo0apTVr1ig9PV0JCQk6d+6cHn/88XLXSUxM1KxZs2paPgAAAAAXsRiGYVR34ezsbMXExGjt2rXq3r27JKlPnz6Kjo7W/Pnzy10nNzdXdrtdycnJuvnmm8tdpkOHDjp9+rQyMzNltVolSc8995z++te/Kjc3t9x1CgoKVFBQ4Hydn5+vsLAw5eXlKTAwsLqHBAAAPCg/P19BQUF8fgP1VI16JlJTU3Xw4EH16NHDOa2wsFAbNmzQSy+9pIKCAmcYKGGz2WS325Wenl7hdm02mxo1alRq3U6dOunHH3/UmTNnyr2A28fHp8KeDgAAAADuV6Mwcf3112vnzp2lpo0ZM0aRkZGaMmVKmSAhSYcOHVJ2drZsNluF2+3du7fefvttFRUVycur+DKOvXv3ymazVflLUAAAAAA8o0YXYAcEBCgqKqrUw9/fXy1atFBUVJSOHz+uyZMna8uWLcrKytL69es1ZMgQtWzZUsOGDXNu584779TUqVOdr++55x4dOnRIEyZM0N69e/XBBx/oqaeeUkJCguuOFAAAAIBLmb7PRHmsVqt27typZcuW6ejRo7LZbOrbt6+WL19e6p4RBw4ccPZASFJYWJjWrl2rSZMmqVu3bmrdurUmTJigKVOmuLI8AAAAAC5Uowuw6zIu4AIAoP7h8xuo30zdZwIAAAAACBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwkQ15eRI69YVP8M9cnJytG7dOuXQyAAAAPUCYaIaHA7JbpdiY4ufHQ5PV9TwOBwO2e12xcbGym63y0Ejuw3B2P0IxgCAXwvCRBVycqRx46SiouLXRUVSfDwnYq6Uk5OjcePGqeh/jVxUVKT4+HhOxNyAYOx+BOPaQzCuHYRjAJXx9nQBdV16+i9BokRhoZSRIYWGeqamhiY9Pd0ZJEoUFhYqIyNDoTSyy1QUjAcM4L3sKhUF4wEDBvBedjGH45f3s5eXlJQkxcV5uqqGx+FwON/TXl5eSkpKUpwHGrqoqEhnzpyp9f0Cv1aNGjWS1Wqt1rKEiSq0b1/8QXX+ua7VKkVEeK6mhqZ9+/by8vIqFSisVqsiaGSXIhi7H8G4dhCMa0ddCcdnzpxRZmZmmb8tAO7VrFkzhYSEyGKxVLocYaIKoaHF33jFxxefeFmt0sKFfGC5UmhoqJKSkhQfH6/CwkJZrVYtXLiQky8XIxi7H8G4dhCMa0ddCMeGYSg3N1dWq1VhYWHy8mJ0NuBuhmHo5MmTOnjwoCTJZrNVujxhohri4oq/8crIKD7x4sPK9eLi4jRgwABlZGQoIiKCIOEGBGP3IxjXDoJx7agL4fjcuXM6efKkLrvsMvn5+dXafoFfuyZNmkiSDh48qODg4EqHPFkMwzBqqzB3ys/PV1BQkPLy8hQYGOjpcoA6KyeHYOxuOTk5BGM3czjKBmOumXA9h8NRJhy7+pqJyj6/T58+rczMTIWHhztPbgDUjlOnTikrK0tt2rSRr69vhcsRJgAA9RLBuHa4OxxXJ0xUdTIDwPWq+/fHMCcAQL0UGkqIqA2hoaH0sAGoEFcyAQAANCB9+vTRxIkTPV1GvWCxWLR69WpPl+FyS5cuVbNmzSpdZubMmYqOjr7ofREmAAAAXKiik/nVq1dX+TObrrBy5UrNmTOn2stnZWXJYrEoLS3NfUXVUbm5uRo0aFCt7W/cuHGyWq1KTk52635GjhypvXv3unUfJQgTAAAADUjz5s0VEBDg6TLqhZCQEPn4+NTKvk6ePKnly5froYceksPhcNt+zp49qyZNmig4ONht+zgfYQIAADR8OTnSunXFz3VUSQ9BcnKyrrnmGvn6+qpLly5av359qeU+//xzXXXVVfLx8ZHNZtMjjzyic+fOOedf2DMSHh6up556SnfddZcCAgJ0+eWXKykpyTm/TZs2kqQrrrhCFotFffr0kSStX79eV111lfz9/dWsWTP17t1b+/fvr7D+KVOmqEOHDvLz81Pbtm01ffp0nT171jl/x44d6tu3rwICAhQYGKgePXroq6++kiTt379fQ4YM0SWXXCJ/f3916dJFa9askVT+kJ0Le3lKhuwsXrxYl19+uZo2bap77rlHhYWFmjt3rkJCQhQcHKwnn3yy1HbOH+ZU0v4rV65U37595efnp+7du2vLli2l1lm0aJHCwsLk5+enYcOG6bnnnqtySJEkrVixQp07d9bUqVO1adMmZWVlVbnO7t279dvf/la+vr7q3Lmz/vWvf5Vb87vvvqs+ffrI19dXb775Zrlt9vTTT+vSSy9VQECA4uLidPr06Sr3Xx2ECQAA0LA5HJLdLsXGFj+78VthV3jooYf04IMPavv27brmmmv0hz/8QYcOHZIkff/99/r973+vK6+8Ujt27NCrr74qh8OhJ554otJtPvvss4qJidH27dt177336p577tHu3bslSVu3bpUk/etf/1Jubq5Wrlypc+fOaejQobruuuv09ddfa8uWLRo3blylw7QCAgK0dOlSffvtt3rhhRe0aNEiPf/88875o0aNUmhoqFJSUpSamqpHHnlEjRo1kiQlJCSooKBAGzZs0M6dO/XMM8+oadOmNWq3ffv26cMPP9RHH32kd955R4sXL9bgwYOVk5Ojzz//XM8884wee+wx/fvf/650O9OmTdPkyZOVlpamDh066LbbbnOGtU2bNmn8+PGaMGGC0tLS1L9//zIBpSIOh0N33HGHgoKC9Pvf/15LliypdPmioiINHTpUfn5++vLLL5WUlKRp06aVu+yUKVN0//33a9euXRowYECZ+e+++65mzJihJ598Ul999ZVsNpteeeWVatVdJaOByMvLMyQZeXl5ni4FAABUU2Wf36dOnTK+/fZb49SpU+Z3kJ1tGF5ehiH98rBai6e7yXXXXWdMmDChzPRVq1YZlZ16ZWZmGpKMp59+2jnt7NmzRmhoqPHMM88YhmEYjz76qNGxY0ejqKjIuczLL79sNG3a1CgsLCx3/3a73bjjjjucr4uKiozg4GDj1VdfLbXf7du3O5c5dOiQIclYv359jY79fHPnzjV69OjhfB0QEGAsXbq03GW7du1qzJw5s9x5S5YsMYKCgkpNu7AtZ8yYYfj5+Rn5+fnOaQMGDDDCw8Od7WIYhtGxY0cjMTHR+VqSsWrVKsMwfmmH//u//3PO/+abbwxJxq5duwzDMIyRI0cagwcPLlXLqFGjytR3ob179xqNGjUy/vvf/zrrDwsLK1XbhT788EPD29vbyM3NdU775JNPyq15/vz5pda9sM169epljB8/vtQyV199tdG9e/cK91/dvz96JgAAQMOVnl76dulS8d0OMzI8U0819OrVy/nf3t7eiomJ0a5duyRJu3btUq9evUr1EPTu3VvHjx9XTiVDuLp16+b8b4vFopCQEB08eLDC5Zs3b67Ro0drwIABGjJkiF544QXl5uZWWvff/vY3/fa3v1VISIiaNm2q6dOn68CBA875DzzwgO6++27169dPTz/9tPbt2+ecd//99+uJJ55Q7969NWPGDH399deV7qs84eHhpa4VufTSS9W5c2d5eXmVmlbZcUul28pms0mSc509e/boqquuKrX8ha/L43A4NGDAALVs2VKS9Pvf/14nTpzQv/71L0nSU089paZNmzofBw4c0J49exQWFqaQkJAq9xUTE1Pp/kveN+e78LVZhAkAANBwtW8veV1wumO1Ft/t0E0CAwOVl5dXZvrRo0dN31i3JDwYhlFmqJHxv/sPVzYEqWQ40fnbK7owZF1gyZIl2rJli6655hotX75cHTp0qHCI0L///W/deuutGjRokN5//31t375d06ZN05kzZ5zLzJw5U998840GDx6szz77TJ07d9aqVaskSXfffbe+++47/elPf9LOnTsVExOjF198UZLk5eXlPMYS51+LUdkxmjnu89cpadOSdSpr/4oUFhZq2bJl+uCDD+Tt7S1vb2/5+fnp8OHDzguxx48fr7S0NOfjsssuK3dfFfH396/Wcu5AmAAAAA1XaKiUlFQcIKTi54UL3XrHw8jISOeFxedLSUlRx44dq1z//BP2c+fOKTU1VZGRkZKkzp07a/PmzaVOYDdv3qyAgAC1bt3aVL2NGzeWVHzSe6ErrrhCU6dO1ebNmxUVFaW333673G1s2rRJdrtd06ZNU0xMjNq3b1/uxdodOnTQpEmTtHbtWt18882lrhsICwvT+PHjtXLlSj344INatGiRJKlVq1Y6duyYTpw44VzWUz9jGxkZ6bzGpER5/6/Pt2bNGh07dkzbt28vFRhWrFih1atX69ChQ2revLkiIiKcD29vb0VGRurAgQP66aefnNtKSUkxVXenTp3KBMGqrh2pLsIEAABo2OLipKys4l9zysoqfu1G9957r/bt26eEhATt2LFDe/fu1csvvyyHw6GHHnqoyvVffvllrVq1Srt371ZCQoKOHDmiu+66y7nt7Oxs3Xfffdq9e7f+/ve/a8aMGXrggQdKDeepieDgYDVp0kQfffSRfvrpJ+Xl5SkzM1NTp07Vli1btH//fq1du1Z79+5Vp06dyt1GRESEDhw4oOTkZO3bt08LFixw9jpI0qlTp/SXv/xF69ev1/79+7Vp0yalpKQ4tzdx4kR9/PHHyszM1LZt2/TZZ58551199dXy8/PTo48+qoyMDL399ttaunSpqWO9WPfdd5/WrFmj5557Tunp6Vq4cKE+/PDDSnsQHA6HBg8erO7duysqKsr5GD58uFq1aqU333yz3PX69++vdu3a6c9//rO+/vprbdq0yXkBdk3vVzJhwgQtXrxYixcv1t69ezVjxgx98803NdpGRQgTAACg4QsNlfr0cWuPRInw8HBt3LhR+/bt0w033KArr7xSS5cu1dKlSzVixIgq13/66af1zDPPqHv37tq4caP+/ve/O8fat27dWmvWrNHWrVvVvXt3jR8/XnFxcXrsscdM1+vt7a0FCxZo4cKFuuyyy3TTTTfJz89Pu3fv1vDhw9WhQweNGzdOf/nLXxQfH1/uNm666SZNmjRJf/nLXxQdHa3Nmzdr+vTpzvlWq1WHDh3SnXfeqQ4dOuiWW27RoEGDNGvWLEnFvSIJCQnq1KmTBg4cqI4dOzp/bah58+Z68803tWbNGnXt2lXvvPOOZs6cafp4L0bv3r312muv6bnnnlP37t310UcfadKkSfL19S13+Z9++kkffPCBhg8fXmaexWLRzTffXOE9J6xWq1avXq3jx4/ryiuv1N133+38/1zR/ioycuRIPf7445oyZYp69Oih/fv365577qnRNipiMaoa6FVP5OfnKygoSHl5eabHIwIAgNpV2ef36dOnlZmZqTZt2tT45Kk+ysrKUps2bbR9+3ZFR0d7uhxU09ixY7V7925t3LjR7fvatGmTfvvb3yojI0Pt2rVz676q+/fn7dYqAAAAgAZk3rx56t+/v/z9/fXhhx/q9ddfd909Gy6watUqNW3aVO3bt1dGRoYmTJig3r17uz1I1ARhAgAAAKimrVu3au7cuTp27Jjatm2rBQsW6O6773bLvo4dO6aHH35Y2dnZatmypfr166dnn33WLfsyi2FOAADAYxjmBNRN1f374wJsAAAAAKYQJgAAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJhCmAAAAGggZs6cyd2zq6lPnz6aOHGip8twuaysLFksFqWlpVW4zPr162WxWHT06NGL3h9hAgAAwEWGDBmifv36lTtvy5Ytslgs2rZtm9v2P3nyZH366ac1Wic8PFzz5893T0F12MqVKzVnzpxa29/bb78tq9Wq8ePHu3U/YWFhys3NVVRUlFv3U4IwAQAA4CJxcXH67LPPtH///jLzFi9erOjoaP3mN79x2/6bNm2qFi1auG37DUnz5s0VEBBQa/tbvHixHn74YSUnJ+vkyZNu2ceZM2dktVoVEhIib29vt+zjQoQJAADQ4OXkSOvWFT+704033qjg4GAtXbq01PSTJ09q+fLliouLq3Dd8PBwzZkzR7fffruaNm2qyy67TC+++GKpZQ4cOKCbbrpJTZs2VWBgoG655Rb99NNPzvkXDnMaPXq0hg4dqnnz5slms6lFixZKSEjQ2bNnJRUP9dm/f78mTZoki8Uii8UiSdq/f7+GDBmiSy65RP7+/urSpYvWrFlTYe1vvvmmYmJiFBAQoJCQEN1+++06ePCgc/6RI0c0atQotWrVSk2aNFH79u21ZMkSScUnwH/5y19ks9nk6+ur8PBwJSYmSip/yM7Ro0dlsVi0fv16Sb8M2fn44491xRVXqEmTJoqNjdXBgwf14YcfqlOnTgoMDNRtt91W6iT+wmFO4eHheuqpp3TXXXcpICBAl19+uZKSkkod5+bNmxUdHS1fX1/FxMRo9erVVQ4pKjmOzZs365FHHlFkZKT+9re/Vbq8JB07dkyjRo2Sv7+/bDabnn/++XJrfuKJJzR69GgFBQVp7Nix5bbZmjVr1KFDBzVp0kR9+/ZVVlZWlfuvLsIEAABo0BwOyW6XYmOLnx0O9+3L29tbd955p5YuXSrDMJzTV6xYoTNnzmjUqFGVrv/Xv/5V3bp107Zt2zR16lRNmjRJn3zyiSTJMAwNHTpUhw8f1ueff65PPvlE+/bt08iRIyvd5rp167Rv3z6tW7dOr7/+upYuXeoMOytXrlRoaKhmz56t3Nxc5ebmSpISEhJUUFCgDRs2aOfOnXrmmWfUtGnTCvdx5swZzZkzRzt27NDq1auVmZmp0aNHO+dPnz5d3377rT788EPt2rVLr776qlq2bClJWrBggf7xj3/o3Xff1Z49e/Tmm28qPDy80mMqz8yZM/XSSy9p8+bNys7O1i233KL58+fr7bff1gcffKBPPvmkTDi70LPPPquYmBht375d9957r+655x7t3r1bUvHJ/ZAhQ9S1a1dt27ZNc+bM0ZQpU6pV2+LFizV48GAFBQXpjjvukKMab8IHHnhAmzZt0j/+8Q998skn2rhxY7lD5P76178qKipKqampmj59epn52dnZuvnmm/X73/9eaWlpuvvuu/XII49Uq+5qMRqIvLw8Q5KRl5fn6VIAAEA1Vfb5ferUKePbb781Tp06ZXr72dmG4eVlGNIvD6u1eLq77Nq1y5BkfPbZZ85p1157rXHbbbdVup7dbjcGDhxYatrIkSONQYMGGYZhGGvXrjWsVqtx4MAB5/xvvvnGkGRs3brVMAzDmDFjhtG9e3fn/D//+c+G3W43zp0755w2YsQIY+TIkaX2+/zzz5fab9euXY2ZM2dW74DLsXXrVkOScezYMcMwDGPIkCHGmDFjyl32vvvuM2JjY42ioqIy8zIzMw1Jxvbt253Tjhw5Ykgy1q1bZxiGYaxbt86QZPzrX/9yLpOYmGhIMvbt2+ecFh8fbwwYMMD5+rrrrjMmTJjgfG2324077rjD+bqoqMgIDg42Xn31VcMwDOPVV181WrRoUer9uGjRojL1XaiwsNAICwszVq9ebRiGYfz3v/81GjVqZKSnp1e4Tn5+vtGoUSNjxYoVzmlHjx41/Pz8ytQ8dOjQUute2GZTp041OnXqVKp9p0yZYkgyjhw5UmEN1f37o2cCAAA0WOnpUlFR6WmFhVJGhvv2GRkZqWuuuUaLFy+WJO3bt08bN27UXXfdVeW6vXr1KvN6165dkqRdu3YpLCxMYWFhzvmdO3dWs2bNnMuUp0uXLrJarc7XNput1BCk8tx///164okn1Lt3b82YMUNff/11pctv375dN910k+x2uwICAtSnTx9JxcOyJOmee+5RcnKyoqOj9fDDD2vz5s3OdUePHq20tDR17NhR999/v9auXVvpvirSrVs3539feuml8vPzU9u2bUtNq+q4z9+GxWJRSEiIc509e/aoW7du8vX1dS5z1VVXVVnX2rVrdeLECQ0aNEiS1LJlS91www3O98fGjRvVtGlT5+Ott97Sd999p7Nnz5baflBQkDp27Fhm+zExMZXuf9euXerZs6dzCJtU9n12MQgTAACgwWrfXvK64GzHapUiIty737i4OL333nvKz8/XkiVLZLfbdf3115vaVslJoGEYpU4IS1Q0vUSjRo3KbK/owoR1gbvvvlvfffed/vSnP2nnzp2KiYmpcIjQiRMndMMNN6hp06Z68803lZKSolWrVkkqHv4kSYMGDdL+/fs1ceJE/fDDD7r++us1efJkSdJvfvMbZWZmas6cOTp16pRuueUW/fGPf5Qkef3vf55x3pCxkus9KjtOi8Vi6rgrW6e8dj6/roosXrxYhw8flp+fn7y9veXt7a01a9bo9ddfV2FhoWJiYpSWluZ8/OEPf3Butzr78/f3r3T/1anxYhAmAABAgxUaKiUlFQcIqfh54cLi6e50yy23yGq16u2339brr7+uMWPGVHrCX+Lf//53mdeRkZGSinshDhw4oOzsbOf8b7/9Vnl5eerUqZPpWhs3bqzCwsIy08PCwjR+/HitXLlSDz74oBYtWlTu+rt379bPP/+sp59+Wr/73e8UGRlZbg9Aq1atNHr0aL355puaP39+qYubAwMDNXLkSC1atEjLly/Xe++9p8OHD6tVq1aS5LyWQ1KVFzu7S2RkpL7++msVFBQ4p3311VeVrnPo0CH9/e9/V3JycqnAkJaWpuPHj+vDDz9UkyZNFBER4XwEBASoXbt2atSokbZu3ercVn5+vtLT02tcd+fOnct9X7lK7fxmFAAAgIfExUkDBhQPbYqIcH+QkIp/onXkyJF69NFHlZeXV+pi5Mps2rRJc+fO1dChQ/XJJ59oxYoV+uCDDyRJ/fr1U7du3TRq1CjNnz9f586d07333qvrrruuyqEulQkPD9eGDRt06623ysfHRy1bttTEiRM1aNAgdejQQUeOHNFnn31WYWC5/PLL1bhxY7344osaP368/vOf/5S5f8Pjjz+uHj16qEuXLiooKND777/v3N7zzz8vm82m6OhoeXl5acWKFQoJCVGzZs3k5eWlnj176umnn1Z4eLh+/vlnPfbYY6aP9WLcfvvtmjZtmsaNG6dHHnlEBw4c0Lx58ySV7UEo8cYbb6hFixYaMWKEs5elxI033iiHw6Ebb7yxzHoBAQH685//rIceekjNmzdXcHCwZsyYIS8vr2qF0vONHz9ezz77rB544AHFx8crNTW1zK+NXYyL6plITEyUxWIp9RNVo0ePdv60WMmjZ8+elW5n6dKlZdaxWCw6ffr0xZQHAAAgqThA9OlTO0GiRFxcnI4cOaJ+/frp8ssvr9Y6Dz74oFJTU3XFFVdozpw5evbZZzVgwABJxSesq1ev1iWXXKJrr71W/fr1U9u2bbV8+fKLqnP27NnKyspSu3btnD0BhYWFSkhIUKdOnTRw4EB17NhRr7zySrnrt2rVSkuXLtWKFSvUuXNnPf30086T7BKNGzfW1KlT1a1bN1177bWyWq1KTk6WVBy8nnnmGcXExOjKK69UVlaW1qxZ4zz5Xrx4sc6ePauYmBhNmDBBTzzxxEUdr1mBgYH65z//qbS0NEVHR2vatGl6/PHHJanUdRTnW7x4sYYNG1YmSEjS8OHD9f7775f6ad/zPffcc+rVq5duvPFG9evXT71791anTp0q3FdFLr/8cr333nv65z//qe7du+u1117TU089VaNtVMZimBxIlZKSoltuuUWBgYHq27ev886Jo0eP1k8//eT87WCp+A3UvHnzCre1dOlSTZgwQXv27Ck1PSQkpNr15OfnKygoSHl5eQoMDKzZwQAAAI+o7PP79OnTyszMVJs2bWp8AlUfhYeHa+LEiaW+pEXd9tZbb2nMmDHKy8tTkyZN3LqvEydOqHXr1nr22WcrvV+Jq1T378/UMKfjx49r1KhRWrRoUbnp0MfHp0ZBQPrlinkAAACgLlq2bJnatm2r1q1ba8eOHZoyZYpuueUWtwSJ7du3a/fu3brqqquUl5en2bNnS5Juuukml+/rYpga5pSQkKDBgwerX79+5c5fv369goOD1aFDB40dO7bKn+GSigOK3W5XaGiobrzxRm3fvt1MaQAAAIBb/Pjjj7rjjjvUqVMnTZo0SSNGjChzl2xXmjdvnrp3765+/frpxIkT2rhxo/Nmf3VFjXsmkpOTtW3bNqWkpJQ7f9CgQRoxYoTsdrsyMzM1ffp0xcbGKjU1VT4+PuWuExkZqaVLl6pr167Kz8/XCy+8oN69e2vHjh1q3759uesUFBSUupo+Pz+/pocCAABQZ2RlZXm6BFTh4Ycf1sMPP1wr+7riiiuUmppaK/u6GDUKE9nZ2ZowYYLWrl1b4dip82/pHhUVpZiYGNntdn3wwQe6+eaby12nZ8+epS7S7t27t37zm9/oxRdf1IIFC8pdJzExUbNmzapJ+QAAAABcqEbDnFJTU3Xw4EH16NHDedONzz//XAsWLJC3t3e5v1Fss9lkt9tr9Lu4Xl5euvLKKytdZ+rUqcrLy3M+zv/NZQAAAADuV6Oeieuvv147d+4sNW3MmDGKjIzUlClTSt2qvcShQ4eUnZ0tm81W7f0YhqG0tDR17dq1wmV8fHwqHDYFAAAAwP1qFCYCAgIUFRVVapq/v79atGihqKgoHT9+XDNnztTw4cNls9mUlZWlRx99VC1bttSwYcOc69x5551q3bq1EhMTJUmzZs1Sz5491b59e+Xn52vBggVKS0vTyy+/7IJDBAAAAOAOLr0DttVq1c6dO7Vs2TIdPXpUNptNffv21fLlyxUQEOBc7sCBA6Vu3nH06FGNGzdOP/74o4KCgnTFFVdow4YNuuqqq1xZHgAAAAAXMn3TurqGm9YBAFD/cNM6oG6q7t+fqftMAAAAoO6YOXOmoqOjK11m9OjRGjp0qNtrsVgsWr16tdv3g7qBMAEAAOBiP/74o+677z61bdtWPj4+CgsL05AhQ/Tpp5+6ZX+TJ09227ZrKjc3V4MGDar28kuXLlWzZs3cVxDcyqXXTAAAAPzaZWVlqXfv3mrWrJnmzp2rbt266ezZs/r444+VkJCg3bt3u2xfhmGosLBQTZs2VdOmTV223YsREhLi6RJQi+iZAAAADV5OTo7WrVunnJwct+/r3nvvlcVi0datW/XHP/5RHTp0UJcuXfTAAw/o3//+d6Xrbt68WdHR0fL19VVMTIxWr14ti8WitLQ0SdL69etlsVj08ccfKyYmRj4+Ptq4cWOZYU6FhYV64IEH1KxZM7Vo0UIPP/ywqrpMtqSHYPXq1erQoYN8fX3Vv3//MvfyevXVV9WuXTs1btxYHTt21BtvvFFq/vnDnLKysmSxWLRy5Ur17dtXfn5+6t69u7Zs2eI8njFjxigvL08Wi0UWi0UzZ86UJL3yyitq3769fH19demll+qPf/xjFS0PTyBMAACABs3hcMhutys2NlZ2u10Oh8Nt+zp8+LA++ugjJSQkyN/fv8z8yobzHDt2TEOGDFHXrl21bds2zZkzR1OmTCl32YcffliJiYnatWuXunXrVmb+s88+q8WLF8vhcOiLL77Q4cOHtWrVqirrP3nypJ588km9/vrr2rRpk/Lz83Xrrbc6569atUoTJkzQgw8+qP/85z+Kj4/XmDFjtG7dukq3O23aNE2ePFlpaWnq0KGDbrvtNp07d07XXHON5s+fr8DAQOXm5io3N1eTJ0/WV199pfvvv1+zZ8/Wnj179NFHH+naa6+tsn7UPoY5AQCABisnJ0fjxo1TUVGRJKmoqEjx8fEaMGCAQkNDXb6/jIwMGYahyMjIGq/71ltvyWKxaNGiRfL19VXnzp31/fffa+zYsWWWnT17tvr371/htubPn6+pU6dq+PDhkqTXXntNH3/8cZU1nD17Vi+99JKuvvpqSdLrr7+uTp06aevWrbrqqqs0b948jR49Wvfee68kOXtb5s2bp759+1a43cmTJ2vw4MGSiu8v1qVLF2VkZCgyMlJBQUGyWCylhkcdOHBA/v7+uvHGGxUQECC73a4rrriiyvpR++iZAAAADVZ6erozSJQoLCxURkaGW/ZXMpTIYrFUutz48eOd1zmUXOuwZ88edevWrdTPcFZ0z62YmJgKt52Xl6fc3Fz16tXLOc3b27vSdSpaLjIyUs2aNdOuXbskSbt27VLv3r1LrdO7d2/n/Iqc33tis9kkSQcPHqxw+f79+8tut6tt27b605/+pLfeeksnT56ssn7UPsIEAABosNq3b1/qRrlS8U12IyIi3LY/i8VS5cn17NmzlZaW5nxIxUHkwhBS0XUO5Q2hcpXygtD508qrsarw1KhRozLrXxjyzhcQEKBt27bpnXfekc1m0+OPP67u3bvr6NGj1TkE1CLCBAAAaLBCQ0OVlJQkq9UqqThILFy40C1DnCSpefPmGjBggF5++WWdOHGizPySk+Hg4GBFREQ4H1JxL8DXX3+tgoIC5/JfffVVjWsICgqSzWYrdbH3uXPnlJqaWuW6586dK7XPPXv26OjRo85hW506ddIXX3xRap3NmzerU6dONa6zROPGjVVYWFhmure3t/r166e5c+fq66+/VlZWlj777DPT+4F7ECYAAECDFhcXp6ysLK1bt05ZWVmKi4tz6/5eeeUVFRYW6qqrrtJ7772n9PR07dq1SwsWLCg19OhCt99+u4qKijRu3Djt2rVLH3/8sebNmyep6mFTF5owYYKefvpprVq1Srt379a9995brW/1GzVqpPvuu09ffvmltm3bpjFjxqhnz57O4VYPPfSQli5dqtdee03p6el67rnntHLlSk2ePLlG9Z0vPDxcx48f16effqqff/5ZJ0+e1Pvvv68FCxYoLS1N+/fv17Jly1RUVKSOHTua3g/cgzABAAAavNDQUPXp08dtPRLna9OmjbZt26a+ffvqwQcfVFRUlPr3769PP/1Ur776aoXrBQYG6p///KfS0tIUHR2tadOm6fHHH5ekUtdRVMeDDz6oO++8U6NHj1avXr0UEBCgYcOGVbmen5+fpkyZottvv129evVSkyZNlJyc7Jw/dOhQvfDCC/rrX/+qLl26aOHChVqyZIn69OlTo/rOd80112j8+PEaOXKkWrVqpblz56pZs2ZauXKlYmNj1alTJ7322mt655131KVLF9P7gXtYjKp+dLieyM/PV1BQkPLy8hQYGOjpcgAAQDVU9vl9+vRpZWZmqk2bNjU+mW4o3nrrLed9GJo0aeLWfS1dulQTJ07kugRIqv7fHz8NCwAAUEcsW7ZMbdu2VevWrbVjxw5NmTJFt9xyi9uDBGAWYQIAAKCO+PHHH/X444/rxx9/lM1m04gRI/Tkk096uiygQgxzAgAAHsMwJ6Buqu7fHxdgo27IyZHWrSt+BgAAQL1AmIDnORyS3S7FxhY/Oxyergi4OIRjwKUayCAKoF6p7t8dYQKelZMjjRsnldwFs6hIio/nJAz1F+EYDYmHg3HJjebOnDnjkf0Dv2YnT56UVPru5eXhAmx4Vnr6L0GiRGGhlJEh1cJvgQMuVVE4HjCA9zPqH4fjl/ezl5eUlCS5+WZvF/L29pafn5/++9//qlGjRvLy4jtQwN0Mw9DJkyd18OBBNWvWzBnqK0KYgGe1b1/8IXV+oLBapYgIz9UEmEU4RkNRR4KxxWKRzWZTZmam9u/fX2v7BSA1a9ZMISEhVS5HmIBnhYYWf9sVH1980mW1SgsXcuKF+olwjIaiDgXjxo0bq3379gx1AmpRo0aNquyRKEGYgOfFxRV/25WRUXzSRZBAfUU4RkNRx4Kxl5cXPw0L1FGECdQNoaGccKFhIByjISAYA6gmbloHAADKl5Pj9mDM5zdQv9EzAQAAykevMYAq8BtrAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEwhTAAAAAAwhTABAAAAwBTCBAAAAABTCBMAAAAATCFMAAAAADCFMAEAAADAFMIEAAAAAFMIEwAAAABMIUwAAAAAMIUwAQAAAMAUwgQAAAAAUwgTAAAAAEy5qDCRmJgoi8WiiRMnOqeNHj1aFoul1KNnz57V3mZycrIsFouGDh16MaUBAAAAcDNvsyumpKQoKSlJ3bp1KzNv4MCBWrJkifN148aNq7XN/fv3a/Lkyfrd735ntiwAAAAAtcRUz8Tx48c1atQoLVq0SJdcckmZ+T4+PgoJCXE+mjdvXuU2CwsLNWrUKM2aNUtt27Y1UxYAAACAWmQqTCQkJGjw4MHq169fufPXr1+v4OBgdejQQWPHjtXBgwer3Obs2bPVqlUrxcXFVauGgoIC5efnl3oAAAAAqD01HuaUnJysbdu2KSUlpdz5gwYN0ogRI2S325WZmanp06crNjZWqamp8vHxKXedTZs2yeFwKC0trdp1JCYmatasWTUtHwAAAICL1ChMZGdna8KECVq7dq18fX3LXWbkyJHO/46KilJMTIzsdrs++OAD3XzzzWWWP3bsmO644w4tWrRILVu2rHYtU6dO1QMPPOB8nZ+fr7CwsBocDQAAAICLUaMwkZqaqoMHD6pHjx7OaYWFhdqwYYNeeuklFRQUyGq1llrHZrPJbrcrPT293G3u27dPWVlZGjJkiHNaUVFRcXHe3tqzZ4/atWtXZj0fH58KezoAAAAAuF+NwsT111+vnTt3lpo2ZswYRUZGasqUKWWChCQdOnRI2dnZstls5W4zMjKyzDYfe+wxHTt2TC+88AK9DQAAAEAdVaMwERAQoKioqFLT/P391aJFC0VFRen48eOaOXOmhg8fLpvNpqysLD366KNq2bKlhg0b5lznzjvvVOvWrZWYmChfX98y22zWrJkklZkOAAAAoO4wfZ+J8litVu3cuVPLli3T0aNHZbPZ1LdvXy1fvlwBAQHO5Q4cOCAvL26+DQAAANRnFsMwDE8X4Qr5+fkKCgpSXl6eAgMDPV0OAACoBj6/gfqN7gEAAAAAphAmAAAAAJhCmAAAAABgCmECAAAAgCmECQAAAACmECYAAAAAmEKYAAAAAGAKYQIAAACAKYSJKuTkSOvWFT8DAAAA+AVhohIOh2S3S7Gxxc8Oh6crarhycnK0bt065ZDa3IZgDAAAXI0wUYGcHGncOKmoqPh1UZEUH8+JmDs4HA7Z7XbFxsbKbrfLQWpzOYJx7SEYux/BGADqDsJEBdLTfwkSJQoLpYwMz9TTUOXk5GjcuHEq+l9jFxUVKT4+nhMxFyIY1x6CsfsRjGsPwRhAdRAmKtC+veR1QetYrVJEhGfqaajS09OdQaJEYWGhMkhtLkMwrh0EY/cjGNcegjGA6iJMVCA0VEpKKg4QUvHzwoXF0+E67du3l9cFqc1qtSqC1OYyBOPaQTB2P4Jx7SAYA6gJwkQl4uKkrKzisblZWcWv4VqhoaFKSkqS9X+pzWq1auHChQoltbkMwbh2EIzdj2BcOwjGAGrCYhiG4ekiXCE/P19BQUHKy8tTYGCgp8tBDeXk5CgjI0MREREECTfJySn+BjcigiDhLg6HQ/Hx8SosLHQG4zi+hXAph6N4aFNh4S/BmCZ2rZycHNnt9lKBwmq1Kisryy3/PvP5DdRvhAkAcCGCsfsRjN2vNoMxn99A/UaYAAAAZdRWMObzG6jfvD1dAAAAqHtCQ0PpXQNQJS7ABgAAAGAKYQIAAACAKYQJAAAAAKYQJgAAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJhCmAAAAABgCmECAAAAgCmECQAAAACmECYAAAAAmEKYAAAAAGAKYQIAAACAKYQJAAAAAKYQJgAAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJhCmAAAAABgCmECAAAAgCmECQAAAACmECYAAAAAmEKYAAAAAGAKYQIAAACAKYQJAAAAAKYQJgAAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJhCmAAAAABgCmECAAAAgCmECQAAAACmECYAAAAAmEKYAAAAAGAKYQIAAACAKYQJAAAAAKYQJgAAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAplxUmEhMTJTFYtHEiROd00aPHi2LxVLq0bNnz0q3s3LlSsXExKhZs2by9/dXdHS03njjjYspDQAAAICbeZtdMSUlRUlJSerWrVuZeQMHDtSSJUucrxs3blzptpo3b65p06YpMjJSjRs31vvvv68xY8YoODhYAwYMMFsiAAAAADcy1TNx/PhxjRo1SosWLdIll1xSZr6Pj49CQkKcj+bNm1e6vT59+mjYsGHq1KmT2rVrpwkTJqhbt2764osvzJQHAAAAoBaYChMJCQkaPHiw+vXrV+789evXKzg4WB06dNDYsWN18ODBam/bMAx9+umn2rNnj6699toKlysoKFB+fn6pBwAAAIDaU+NhTsnJydq2bZtSUlLKnT9o0CCNGDFCdrtdmZmZmj59umJjY5WamiofH58Kt5uXl6fWrVuroKBAVqtVr7zyivr371/h8omJiZo1a1ZNywcAAADgIhbDMIzqLpydna2YmBitXbtW3bt3l1Q8RCk6Olrz588vd53c3FzZ7XYlJyfr5ptvrnDbRUVF+u6773T8+HF9+umnmjNnjlavXq0+ffqUu3xBQYEKCgqcr/Pz8xUWFqa8vDwFBgZW95AAAIAH5efnKygoiM9voJ6qUc9EamqqDh48qB49ejinFRYWasOGDXrppZecvQrns9lsstvtSk9Pr3TbXl5eioiIkCRFR0dr165dSkxMrDBM+Pj4VNrTAQAAAMC9ahQmrr/+eu3cubPUtDFjxigyMlJTpkwpEyQk6dChQ8rOzpbNZqtRYYZhlOp5AAAAAFC31ChMBAQEKCoqqtQ0f39/tWjRQlFRUTp+/Lhmzpyp4cOHy2azKSsrS48++qhatmypYcOGOde588471bp1ayUmJkoqvv4hJiZG7dq105kzZ7RmzRotW7ZMr776qgsOEQAAAIA7mL7PRHmsVqt27typZcuW6ejRo7LZbOrbt6+WL1+ugIAA53IHDhyQl9cvPyR14sQJ3XvvvcrJyVGTJk0UGRmpN998UyNHjnRleQAAAABcqEYXYNdlXMAFAED9w+c3UL+Zus8EAAAAABAmAAAAAJhCmIDn5eRI69YVPwMAAKDeIEzAsxwOyW6XYmOLnx0OT1cEmEcwBgD8yhAm4Dk5OdK4cVJRUfHroiIpPp4TMdRPBGM0JARjANVEmIDnpKf/EiRKFBZKGRmeqQcwi2CMhoRgDKAGCBPwnPbtJa8L3oJWqxQR4Zl6ALMIxmgoCMYAaogwAc8JDZWSkooDhFT8vHBh8XSgPiEYo6EgGAOoIcIEPCsuTsrKKh6bm5VV/BqobwjGaCgIxgBqyNvTBQAKDeWkC/VfXJw0YEDxN7gREbynUT+VBOP4+OIeCYIxgCpYDMMwPF2EK+Tn5ysoKEh5eXkKDAz0dDkAANRfOTm1Foz5/AbqN3omAABAafQYA6gmrpkAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJhCmAAAAABgCmECAAAAgCmECQAAAACmECYAAAAAmEKYAAAAAGAKYQIAAACAKYQJAAAAAKZ4e7oAVzEMQ5KUn5/v4UoAAEB1lXxul3yOA6hfGkyYOHbsmCQpLCzMw5UAAICaOnbsmIKCgjxdBoAashgN5KuAoqIi/fDDDwoICJDFYqly+fz8fIWFhSk7O1uBgYG1UGH9QLuUj3YpH+1SPtqlYrRN+X7N7WIYho4dO6bLLrtMXl6MvgbqmwbTM+Hl5aXQ0NAarxcYGPir+4e7OmiX8tEu5aNdyke7VIy2Kd+vtV3okQDqL74CAAAAAGAKYQIAAACAKb/aMOHj46MZM2bIx8fH06XUKbRL+WiX8tEu5aNdKkbblI92AVBfNZgLsAEAAADUrl9tzwQAAACAi0OYAAAAAGAKYQIAAACAKYQJAAAAAKY0yDCRlZWluLg4tWnTRk2aNFG7du00Y8YMnTlzptzlDx06pNDQUFksFh09erTS7VoslnIfK1ascNPRuI672qXEli1bFBsbK39/fzVr1kx9+vTRqVOnXHwUrufOdunTp0+Z98qtt97qhqNwPXe/X6TiO98OGjRIFotFq1evdl3xbubOtomPj1e7du3UpEkTtWrVSjfddJN2797thqNwPXe1y+HDh3XfffepY8eO8vPz0+WXX677779feXl5bjoS13Ln+yUpKUl9+vRRYGBgjf72AMBVGswdsM+3e/duFRUVaeHChYqIiNB//vMfjR07VidOnNC8efPKLB8XF6du3brp+++/r3S7YWFhys3NLTUtKSlJc+fO1aBBg1x6DO7grnaRioPEwIEDNXXqVL344otq3LixduzYIS+vup9X3dkukjR27FjNnj3b+bpJkyYuq92d3N0ukjR//nxZLBZXll0r3Nk2PXr00KhRo3T55Zfr8OHDmjlzpm644QZlZmbKarW643Bcxl3t8sMPP+iHH37QvHnz1LlzZ+3fv1/jx4/XDz/8oL/97W/uOhyXcef75eTJkxo4cKDz318AqHXGr8TcuXONNm3alJn+yiuvGNddd53x6aefGpKMI0eO1Gi70dHRxl133eWiKmufq9rl6quvNh577DE3VVn7XNUu1113nTFhwgT3FOkBrvw7SktLM0JDQ43c3FxDkrFq1SrXF1yL3PVvzI4dOwxJRkZGhosqrV3uapd3333XaNy4sXH27FkXVVq7XN0u69atM9WOAHCx6v7Xxi6Sl5en5s2bl5r27bffavbs2Vq2bJmpb9BTU1OVlpamuLg4V5VZ61zRLgcPHtSXX36p4OBgXXPNNbr00kt13XXX6YsvvnBX2W7nyvfLW2+9pZYtW6pLly6aPHmyjh075upya42r2uXkyZO67bbb9NJLLykkJMQdpdY6d/wbc+LECS1ZskRt2rRRWFiYq0qtVe5ol5LtBgYGytu7fnawu6tdAKDWeTrN1IaMjAwjMDDQWLRokXPa6dOnjW7duhlvvPGGYRjmvtW55557jE6dOrm63FrjqnbZsmWLIclo3ry5sXjxYmPbtm3GxIkTjcaNGxt79+5192G4nCvfL0lJScYnn3xi7Ny503jnnXeM8PBwo1+/fu4s321c2S7jxo0z4uLinK9Vz3smXP1vzMsvv2z4+/sbkozIyMh62yvhrn97f/75Z+Pyyy83pk2b5uqSa4U72oWeCQCeUq/CxIwZMwxJlT5SUlJKrfP9998bERERpU5cDMMwJk2aZIwcOdL5uqb/EJ88edIICgoy5s2bd9HHdbE83S6bNm0yJBlTp04tNb1r167GI488cvEHaJKn26U8X331lSHJSE1NNX1cF8vT7fL3v//diIiIMI4dO+acVlfChKfbpsTRo0eNvXv3Gp9//rkxZMgQ4ze/+Y1x6tQplxyjGXWlXQzDMPLy8oyrr77aGDhwoHHmzJmLPraLUZfahTABwFMshmEYZno0POHnn3/Wzz//XOky4eHh8vX1lVR80V7fvn119dVXa+nSpaW6jaOjo7Vz507nxZ+GYaioqEhWq1XTpk3TrFmzKt3PG2+8obi4OH3//fdq1arVRR7ZxfF0u2RmZqpt27Z64403dMcddzinjxw5Ut7e3nrrrbdccZg15ul2KY9hGPLx8dEbb7yhkSNHmjyyi+Ppdpk4caIWLFhQajuFhYXy8vLS7373O61fv94FR2mOp9umPGfOnNEll1yi//u//9Ntt91m8sguTl1pl2PHjmnAgAHy8/PT+++/79yfp9SVdpGk9evXq2/fvjpy5IiaNWt2cQcGADVQrwabtmzZUi1btqzWst9//7369u2rHj16aMmSJWXGn7733nulfrY0JSVFd911lzZu3Kh27dpVuX2Hw6E//OEPHg8SkufbJTw8XJdddpn27NlTavrevXs9+itXnm6X8nzzzTc6e/asbDZbtddxNU+3yyOPPKK777671LSuXbvq+eef15AhQ2p4NK7l6bapiGEYKigoqNE6rlQX2iU/P18DBgyQj4+P/vGPf3g8SEh1o10AwOM81CPiViXdyLGxsUZOTo6Rm5vrfFSkvC7inJwco2PHjsaXX35Zatn09HTDYrEYH374obsOwS3c2S7PP/+8ERgYaKxYscJIT083HnvsMcPX17dejPV2V7tkZGQYs2bNMlJSUozMzEzjgw8+MCIjI40rrrjCOHfunLsP66K5++/ofKojw5yqy11ts2/fPuOpp54yvvrqK2P//v3G5s2bjZtuuslo3ry58dNPP7n7sC6au9olPz/fuPrqq42uXbsaGRkZpbb7a/9bys3NNbZv324sWrTIkGRs2LDB2L59u3Ho0CF3HhIAONWrnonqWrt2rTIyMpSRkaHQ0NBS84wajOo6e/as9uzZo5MnT5aavnjxYrVu3Vo33HCDS+qtLe5sl4kTJ+r06dOaNGmSDh8+rO7du+uTTz6pF9+ouatdGjdurE8//VQvvPCCjh8/rrCwMA0ePFgzZsyo8/cLkNz/d1SfuattfH19tXHjRs2fP19HjhzRpZdeqmuvvVabN29WcHCwS4/BHdzVLqmpqfryyy8lSREREaWWzczMVHh4+MUV7mbu/Ft67bXXSg2BuvbaayVJS5Ys0ejRoy+ucACohnp1zQQAAACAuoMfsgYAAABgCmECAAAAgCmECQAAAACmECYAAAAAmEKYAAAAAGAKYQIAAACAKYQJAAAAAKYQJgAAAACYQpgAAAAAYAphAgAAAIAphAkAAAAAphAmAAAAAJjy/wgqgH+ZxaM9AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D # noqa\n", "\n", "fig, ax = plt.subplots()\n", "nind = 3\n", "ax1 = ax.plot(\n", " fieldsetA.U.grid.lon[:nind, :nind],\n", " fieldsetA.U.grid.lat[:nind, :nind],\n", " \".r\",\n", " label=\"U points assuming A-grid\",\n", ")\n", "ax2 = ax.plot(\n", " fieldsetA.V.grid.lon[:nind, :nind],\n", " fieldsetA.V.grid.lat[:nind, :nind],\n", " \".b\",\n", " label=\"V points assuming A-grid\",\n", ")\n", "\n", "ax3 = ax.plot(\n", " fieldsetC.U.grid.lon[:nind, :nind],\n", " fieldsetC.U.grid.lat[:nind, :nind],\n", " \".k\",\n", " label=\"C-grid points\",\n", ")\n", "\n", "ax.legend(handles=[ax1[0], ax2[0], ax3[0]], loc=\"center left\", bbox_to_anchor=(1, 0.5))\n", "plt.show()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Further information about the NEMO C-grids is available in [the NEMO 3D tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_nemo_3D.html).\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Arakawa B-grids\n", "\n", "Interpolation for Arakawa B-grids is detailed in [Section 2.1.3 of Delandmeter and Van Sebille (2019)](https://www.geosci-model-dev.net/12/3571/2019/gmd-12-3571-2019.html). Again, the dimensions that need to be provided are those of the barycentric cell edges `(i, j, k)`.\n", "\n", "To load B-grid data, you can use the method `FieldSet.from_b_grid_dataset`, or specifically in the case of POP-model data `FieldSet.from_pop`.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 3D C- and B-grids\n", "\n", "For 3D C-grids and B-grids, the idea is the same. It is important to follow the indexing notation, which is defined in Parcels and in [Delandmeter and Van Sebille (2019)](https://www.geosci-model-dev.net/12/3571/2019/gmd-12-3571-2019.html). In 3D C-grids, the vertical (`W`) velocities are at the top and bottom. Note that in the vertical, the bottom velocity is often omitted, since a no-flux boundary conditions implies a zero vertical velocity at the ocean floor. That means that the vertical dimension of `W` often corresponds to the amount of vertical levels.\n" ] } ], "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.11.6" } }, "nbformat": 4, "nbformat_minor": 4 }