{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Documentation on dimensions and indexing in Parcels"
   ]
  },
  {
   "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: **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`.** 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](https://oceanparcels.org/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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parcels import FieldSet\n",
    "from glob import glob\n",
    "import numpy as np\n",
    "from os import path\n",
    "from datetime import timedelta as delta"
   ]
  },
  {
   "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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 https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-551c32df110d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m               'W': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthw', 'time': 'time_counter'}}\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mfieldset\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFieldSet\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_nemo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilenames\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvariables\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdimensions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/Codes/PARCELScode/parcels/fieldset.py\u001b[0m in \u001b[0;36mfrom_nemo\u001b[0;34m(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, field_chunksize, **kwargs)\u001b[0m\n\u001b[1;32m    463\u001b[0m         fieldset = cls.from_c_grid_dataset(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,\n\u001b[1;32m    464\u001b[0m                                            \u001b[0mallow_time_extrapolation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mallow_time_extrapolation\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtracer_interp_method\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtracer_interp_method\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 465\u001b[0;31m                                            field_chunksize=field_chunksize, **kwargs)\n\u001b[0m\u001b[1;32m    466\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfieldset\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'W'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    467\u001b[0m             \u001b[0mfieldset\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_scaling_factor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1.\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Codes/PARCELScode/parcels/fieldset.py\u001b[0m in \u001b[0;36mfrom_c_grid_dataset\u001b[0;34m(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, field_chunksize, **kwargs)\u001b[0m\n\u001b[1;32m    524\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    525\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;34m'U'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdimensions\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;34m'V'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdimensions\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mdimensions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'U'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mdimensions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'V'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 526\u001b[0;31m             raise ValueError(\"On a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U and V. \"\n\u001b[0m\u001b[1;32m    527\u001b[0m                              \"See also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb\")\n\u001b[1;32m    528\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;34m'U'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdimensions\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;34m'W'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdimensions\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mdimensions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'U'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mdimensions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'W'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\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 https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb"
     ]
    }
   ],
   "source": [
    "data_path = 'NemoNorthSeaORCA025-N006_data/'\n",
    "ufiles = sorted(glob(data_path+'ORCA*U.nc'))\n",
    "vfiles = sorted(glob(data_path+'ORCA*V.nc'))\n",
    "wfiles = sorted(glob(data_path+'ORCA*W.nc'))\n",
    "\n",
    "filenames = {'U': ufiles, 'V': vfiles, 'W': wfiles}\n",
    "variables = {'U': 'uo', 'V': 'vo', 'W': 'wo'}\n",
    "dimensions = {'U': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthu', 'time': 'time_counter'},\n",
    "              'V': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthv', 'time': 'time_counter'},\n",
    "              'W': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthw', 'time': 'time_counter'}}\n",
    "\n",
    "fieldset = FieldSet.from_nemo(filenames, variables, dimensions)"
   ]
  },
  {
   "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**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fieldsetA = FieldSet.from_netcdf(filenames, variables, dimensions)"
   ]
  },
  {
   "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: File NemoNorthSeaORCA025-N006_data/coordinates.nc could not be decoded properly by xarray (version 0.15.1).\n",
      "         It will be opened with no decoding. Filling values might be wrongly parsed.\n",
      "WARNING: Casting lon data to np.float32\n",
      "WARNING: Casting lat data to np.float32\n",
      "WARNING: Trying to initialize a shared grid with different chunking sizes - action prohibited. Replacing requested field_chunksize with grid's master chunksize.\n"
     ]
    }
   ],
   "source": [
    "mesh_mask = data_path + 'coordinates.nc'\n",
    "\n",
    "filenames = {'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",
    "dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'},\n",
    "              'V': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'},\n",
    "              'W': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'}}\n",
    "\n",
    "fieldsetC = FieldSet.from_nemo(filenames, variables, dimensions)"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiEAAAD4CAYAAAA6uTZJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de1iUZf4/8PdnZhgE8QCCBw6GqcMwgGiy2LhWs8L6tRb9Zmrr4iHX9Zv+zP128FBua+2WLpV15YVdXuW2ZWtmXmVlS9/aVWrKcjbDMyKgbrjmEU8EnoCZ+/fHgIs4CDjg48y8X9c11wPz3M89n5lbmffczzPPI0opEBEREd1oOq0LICIiosDEEEJERESaYAghIiIiTTCEEBERkSYYQoiIiEgTBq0LaCwyMlLFx8drXQYRkU/ZunXrSaVUlNZ1ELXGTRdC4uPjUVBQoHUZREQ+RUQOal0DUWtxdwwRERFpgiGEiIiINMEQQkRERJpgCCEiIiJNMIQQERGRJhhCiIiISBN+F0IcDiAnx70kIiKim9dNd54QbzgcQEYGUF0NGI1Afj5gtWpdVftzOByw2+2w2WywBsATdjgAux2w2QJjfIHAG2MgMMeZKND4VQix290BxOl0L+12///j5XA4kJGRgerqahiNRuTn5/v1m1QgBs1AG2MgcMc50IImkV/tjrHZ3H+w9Hr30mbTuqL2Z7fbUV1dDafTierqatjtdq1Laleegqa/C7QxBgJvnOuD5sKFC5GRkQEH9ydTgPCrEGK1uj8xPftsYHxyAgCbzQaj0Qi9Xg+j0QibnyevQAyagTbGQOCNcyAGTSIAEKWU1jVcIS0tTfHaMa0TaNO4gXisQKCNMRBY49wWu9xEZKtSKq2dSiRqFwwhREQ3AW+DJkMI+SK/OjCViMhXWa3WgJnlIqrnV8eEEBERke9gCCEiIiJNMIQQERGRJhhCiIiISBMMIURERKQJhhAiIiLSBEMIERERaYIhhIiIiDTBEEJERESaaHEIERG9iGwXkby631eKyPcisqPuNtDDNj9rsH6HiFwUkXvb8gkQERGRb2rNadsfBrAXQOcG981TSr3f1AZKqS8ADAQAEYkAsB/AP66jTiIiIvIzLZoJEZFYAL8A8LoXjzUOwKdKqfNe9EFERER+oqW7Y5YCmA/A1ej+xSKyS0ReFpHgZvqYAGCNpxUi8qCIFIhIQXl5eQtLIiIiIl/WbAgRkSwAJ5RSWxutWgDADOAnACIAPH6NPnoBSAHwd0/rlVIrlFJpSqm0qKioltZOREREPqwlMyE/BTBaRMoAvAtguIi8rZQ6qtwuAXgTQPo1+rgfwIdKqRqvKyYiIiK/0GwIUUotUErFKqXi4d6l8rlSalLd7AZERADcC6DwGt38Ck3siiEiIqLA5M15QlaLyG4AuwFEAlgEACKSJiKXD2AVkXgAcQC+9OKxiIiIyM+05iu6UErZAdjrfh7eRJsCANMb/F4GIOZ6CyQiIiL/xDOmEhERkSYYQoiIiEgTDCFERESkCYYQIiIi0gRDCBEREWmCIYSIiIg0wRBCREREmmAIISIiIk0whBAREZEmGEKIiIhIEwwhREREpAmGECIiItIEQwgRERFpgiGEiIiINMEQ4uscDiAnx70k/8QxJiI/ZdC6APKCwwFkZADV1YDRCOTnA1ar1lVRW+IYBwaHA7DbAZuN40sBhTMhvsxud785OZ3upd2udUXU1jjG/q8+aC5c6F5yxosCCEOIL7PZ3J+O9Xr30mbTuiJqaxxj/8egSQGMu2N8mdXqnp7nNK7/4hj7v/qgWb/LjUGTAogopbSu4QppaWmqoKBA6zKIiG6cNjgmRES2KqXS2rQuonbGmRAiIq1ZrZzlooDEY0KIiIhIEwwhREREpAmGECIiItIEQwgRERFpgiGEiIiINMEQQkRERJpgCCEiIiJNMIQQERGRJhhCiIiISBMMIURERKQJhhAiIiLSBEMIERERaYIhhIiIiDTR4hAiInoR2S4ieXW/rxSR70VkR91tYBPb9RaRf4jIXhEpEpH4timdiIiIfJmhFW0fBrAXQOcG981TSr3fzHZ/BbBYKbVBRMIAuFpZIxEREfmhFs2EiEgsgF8AeL01nYuIBYBBKbUBAJRSVUqp862ukoiIiPxOS3fHLAUwH1fPYiwWkV0i8rKIBHvYzgTgrIh8ULcrZ4mI6Bs3EpEHRaRARArKy8tb9wyIiIjIJzUbQkQkC8AJpdTWRqsWADAD+AmACACPe9jcAOAOAHPr2t0KYGrjRkqpFUqpNKVUWlRUVKueABEREfmmlsyE/BTAaBEpA/AugOEi8rZS6qhyuwTgTQDpHrb9AcB2pdS/lFK1AD4CcFsb1U5EREQ+rNkQopRaoJSKVUrFA5gA4HOl1CQR6QUAIiIA7gVQ6GHz7wCEi0j99MZwAEVtUjkRERH5NG/OE7JaRHYD2A0gEsAiABCRNBF5HQCUUk64d8Xk17UVAH/2rmQiIiLyB6KU0rqGK6SlpamCggKtyyAi8ikislUplaZ1HUStwTOmEhERkSYYQoiIiEgTDCFERESkCYYQIiIi0gRDCBEREWmCIYSIiIg0wRBCREREmmAIISIiIk0YtC6AiIjax9atW7sbDIbXASSDHzrpxnMBKKytrZ0+ePDgE54aMIQQEfkpg8Hwes+ePROjoqLO6HS6m+v02OT3XC6XlJeXW44dO/Y6gNGe2vhdMnY4gJwc9zIQOBwO5OTkwBEoTxgcY6JWSI6KivqRAYS0oNPpVFRUVAXcM3Ee+dVMiMMBZGQA1dWA0Qjk5wNWq9ZVtR+Hw4GMjAxUV1fDaDQiPz8fVn9+wuAYB8oY2+2AzebfY9uQw+GA3W6HzWZr6/HVMYCQlur+/TU54eFXMyF2u/vNyel0L+12rStqX3a7HdXV1XA6naiurobd358wOMb+Psb1IXPhQvcyECZ/6oPmwoULkZGRwRkvCih+FUJsNvenY73evbTZtK6ofdlsNhiNRuj1ehiNRtj8/QmDY+zvYxxoIRPw76BZUlJi7N+/f1LD+x577LHop556qoe3fX/11VehU6dOjbtWm5MnT+qfe+65KG8f60YqKysLGjly5K1t2WdGRkbfgQMHmtuir9WrV3f53e9+19PTutDQ0EGt7c+vdsdYre7p+UCZyrVarcjPz2+vadybEsfYv59wfcis393m55kLwH+CZv0uN38Pmm3lzjvvPH/nnXeev1abU6dO6f/yl790f+KJJ8pvVF3eio+Pr/nss8/+1Vb9nTx5Ur9nz56OoaGhzuLiYqPZbK6+3r5qamowceLECgAVbVWfX82EAO43pQUL/P/NqZ7VasWCBQv8/s2pIY6x/6oPmc8+6//H+9SrD5rPPvvszXHMz8aNHbFgQU9s3NjxRj7s2LFj47Ozs3sPHjw4IT4+PnnNmjVdAOD8+fMybty4eJPJZElMTLT87W9/6wQAeXl5nX72s5/1A9yzK+PHj49PT09PiI2NTVm0aFF3AJgzZ07soUOHgs1ms2XGjBmxBw8eDEpLS0swm82W/v37J3322WdhjeuYO3dur+Tk5MT+/fsn/epXv7rF5XIBABYtWtS9b9++SSaTyZKVlXUrAHzyySdhZrPZYjabLYmJiZYzZ87oGtYFAFOmTOmdm5vbDQBiYmJSZs+eHTNw4EBzcnJy4tdffx06bNiw/nFxcckvvPBCFHDl7FFubm63ESNG9L3jjjv633LLLckzZ86Mre/35ZdfjoyPj09OT09PmDBhwi1Tpkzp7el1XbVqVXhmZubZMWPGnH7rrbcimnr9m+pv7Nix8dOnT48dMmSIadasWbG5ubnd6tcVFxcb65/Lww8/HN3iwW7Ar2ZCiMj3Wa2BET4aslqt2ocPwB1AsrJMqKnR4eWXXcjLK0Vm5rkb9fCHDh0K3rJlS0lRUVFwZmZmwn//93/vfv7557sDQGlpadH27ds73HPPPf0PHDhQ2Hjb/fv3d9i8eXPJ2bNn9YmJicnz5s0rf+mll37IysoKKS4uLgKAp59+ukdGRkbF888/f6y2thaVlZVXfRCfN2/eiRdffPEoANx777193n333S7Z2dkVubm5PQ8ePLg7JCREnTx5Ug8AL730Us/c3NyDI0aMOFdRUaELDQ11Nfcc4+Liqnfs2FH8m9/8Jm7atGnx3377bfGFCxd0ycnJSfPnz79qxqaoqCh0586dRSEhIa5+/folz50797jBYMCLL77Ya9u2bUVdu3Z1DR061JSUlHTB0+O99957EU899dSR6OjomnHjxvXNyck51rhNWVlZ0LX6O3DgQIdvvvmm1GAwoD5QAcCsWbN6T58+vXz27NmncnJyrmu3l9/NhBAR0XXKz++EmhodXC6gtlaH/PxO3nQnIq26f+zYsaf1ej1SUlIuxcXFXdqxY0eHzZs3h02ZMuUUAAwaNOhidHR09e7duzs03nbEiBFnQ0JCVK9evWojIiJqfvjhh6s+ZN9+++3n1qxZE/nYY49Fb9myJSQ8PPyq0PDpp592GjBggNlkMlk2b97cqbCwMAQAEhISLowZM6bP8uXLI4KCglRdf1Vz586NW7RoUfeTJ0/qg4KCmn1N7r///rMAkJKScv622247Fx4e7oqOjq4NDg521YebhoYNG/Zjt27dnKGhoapfv34XDxw4ELxp06aOQ4YMqezRo4czODhYjRkz5oynxzp06JDh4MGDwSNGjKgaMGDAJYPBoL777rurXrvm+rvvvvvOGAxXz1ls27Yt7H/+539OA8CMGTNONfvkPWAIISIit4yMSgQFuaDXAwaDCxkZld5016NHj9qKioor3lhPnz6tj4yMrPXUvnE4EREo1bJvGAcHB19uqNfrUVtbe1XSufvuu6u++uqrkpiYmOqpU6f2eeWVV7o1XH/+/HmZM2fOLR988MGB0tLSokmTJp28ePGiDgC++OKLfQ899FD51q1bO6amplpqamrwpz/96djrr79+8MKFC7qhQ4cmbt++vUNQUJCq34UDAJcuXbqijg4dOigA0Ol0MBqNl2vW6XSoqam5quaGbfR6vaqpqZGWviZvvfVWxI8//qiPi4tLiYmJSTl8+HDwqlWrImpra1G/G+mRRx6Jbq6/sLCwJmd4vP0KOEMIERG5ZWaeQ15eKebNO9wWu2K6dOni6t69e8369es7AcDx48f1dru9y/Dhw6s8tf/ggw/CnU4n9uzZE3zo0KHg1NTUi8OGDat6++23IwBg165dwUePHjUOGDDgYgsf33nu3LnL73OlpaXGmJiYmjlz5pycNGnSyW3btoU2bH/+/HkdAPTs2bO2oqJC97e//S0cAJxOJw4cOGAcNWpU5fLly3+orKzUV1RU6Pfs2ROcnp5+YfHixcdSUlLOFRYWdujbt++l/fv3h1y4cEFOnTql//rrrztf36vXtDvuuOPct99+26m8vFxfU1OD9evXh3tq9/7770d8+OGH+w4fPrz78OHDu7/99tuijz76KMJgMKC4uLiouLi4aOnSpUda2l9jt912W9Wf//znCAD485//3K259p7wmBAiIvqPzMxzbXkcyFtvvfX9rFmzej/++ONxAPD4448fSUpKuuSpbb9+/S6lp6cnnDp1Kmjp0qUHQ0ND1fz5809Mnjz5FpPJZNHr9XjttdfKQkJCWvTpu2fPns7BgwdX9e/fP2n48OEVycnJF3Jzc3saDAYVGhrqXL169fcN20dGRjonTpxYbrFYkmJjY6tTU1PPAUBtba1kZ2f3qays1CulZMaMGccjIyOdc+bMid68eXNnnU6nTCbThXHjxlWEhISoUaNGnUlMTEzq06fPxaSkpGt+g+d69OnTp+bRRx89+pOf/CSxe/fuNSaT6UKXLl2cDduUlJQYjxw5Yhw+fPjlsTSbzdVhYWHOzz//vGPD+1vSnyfLly//94QJE25dvnx5j9GjR3vcJdScFk/r3ChpaWmqoKBA6zKIiHyKiGxVSqU1vG/nzp1lqampJ7WqqTXGjh0bn5WVVfHrX//6ut7MAk1FRYWuS5curpqaGvzXf/1Xv6lTp56cMmXK2Zulv4Z27twZmZqaGu9pHXfHEBER+Zh58+ZFm81mi8lkSurdu/elSZMmeRUY2rq/luLuGCIi0ty6devKtK7Bl6xYseKHm7m/luJMCBEREWmCIYSIiIg0wRBCREREmmAIISIiIk0whBARUbtIT09PWLdu3RUn63rmmWe6T5o0yePF1lrjl7/85S1bt2696hTkDa1ataprc21uNo888kj0Rx995NXp8hv65ptvQkRkcONxuF533XVXP0+nl3/sscein3rqqR6t7Y8hhIiI2sX48eNPrVmz5oort65bty5i0qRJp73te+3atQcHDx58zTOnfvTRR1137doV4u1j3UhLly49cu+993p1uvyGVq1a1e22226reuedd5q8gm5LuFwuOJ1OfPnll/sjIyObPZFZSzGEEBHRZRs3ouOCBei5cSM6etvX5MmTz+Tn53e5cOGCAO6zeJ44cSJoxIgRV5y2vaSkxNinT5+k++67L95kMllGjhx5a/0VbtevX98pMTHRYjKZLOPHj4+v7ys9PT3hq6++CgWA0NDQQb/97W9jEhISLKmpqeZDhw4ZNmzY0HHjxo1df//738eazWbLnj17ghctWtS9b9++SSaTyZKVlXVr43pLSkqMgwcPTrBYLIkWiyVxw4YNHQHg4MGDQWlpaQlms9nSv3//pM8++yystrYWY8eOje/fv3+SyWSy/PGPf+zeuK6jR48aYmJiUgAgNze3W2ZmZt/hw4f3i4mJSfnTn/4U9Yc//KFHYmKiJTU11Xz8+HE94D5p25tvvhkOADExMSmPPvpotMViSTSZTJbt27d3AIAjR44Yhg4d2t9isSRmZ2ffEh0dnXL06NGrTrnhcrmQl5cX/te//rVs06ZNnc+fP+/xyoFN9VdSUmK89dZbkyZNmtQ7KSnJcuDAAWNMTMzlx3r88cd7xsfHJw8dOtS0b9++4Nb/C2EIISKiOhs3omNWFkwvvICYrCyYvA0iPXv2dKampp5bt25dF8B9QbXRo0ef0emufuspKyvrMHPmzPLS0tKiTp06uZYsWRJ1/vx5mTFjRp+1a9ceKC0tLaqtrcWSJUuuumT8hQsXdFartaqkpKTIarVWLVu2LOrnP//5uczMzLOLFi36obi4uCgpKelSbm5uz8LCwqLS0tKilStXHmzcT3R0dO2mTZtKi4qK9q5du/Zfjz76aG8AeOONNyIyMjIqiouLi/bu3btnyJAh5x0OR+jRo0eD9u3bt6e0tLTooYceavYqsqWlpSHr1q3713fffbc3JycnJjQ01LV3796itLS0c6+99prHa69ERkbWFhUV7Z02bVr5c8891wMAnnjiiei77rqrsqioaO9999135ujRo0ZP227YsCEsLi7uUlJS0qUhQ4ZUvvfee108tbtWf2VlZR1+/etfn9q7d2+RyWSqrr9/06ZNoR9++GHE7t27i/Ly8vbv3Lnzuv6tMIQQEREAID8fnWpqoHO5gNpa6PLz4fWxCffff//ptWvXhgPABx98EDF58mSPu2J69uxZPWLEiHMAMHny5FObN28O27lzZ4fY2NhLAwYMuAQAU6dOPfX1119fVVNQUJCaMGFCBQAMHjz43MGDBz2+KSckJFwYM2ZMn+XLl0cEBQVddc2S6upqyc7Ojq+bdel74MCBDgBw++23n1uzZk3kY489Fr1ly5aQ8PBwl9lsvnTo0KHgBx54IO7999/vHB4e3uwuiqFDh1aGh4e7oqOja8PCwpzjx48/CwApKSnny8rKPM4kZGdnnwGA9PT084cOHQoGgC1btoQ98MADpwFg3LhxP3bu3NnjY7/99tsR48aNOw0AEyZMOP3uu+963CVzrf569epVnZGRcdW1hL744ouwe+6552ynTp1cERERrhEjRlzXGVYZQoiICACQkYHKoCC49HrAYIArIwNeH5swceLEs998803nr7/+OvTixYu6YcOGebygm4hc9XtLr21mMBhU/eyKwWBAbW2tx90OX3zxxb6HHnqofOvWrR1TU1MtNTU1V6xfvHhxj+7du9fs3bu3aPfu3UU1NTU6ALj77rurvvrqq5KYmJjqqVOn9nnllVe6RUVFOQsLC4t+9rOfVS5fvrz7hAkT4utrcTrd7+GNd38YjcbLT0in06FDhw6q/uemaq5vYzAYVH2blrwutbW1+PTTT8OXLFkSHRMTkzJv3rzeX375ZZczZ87ocnJyosxms8VsNlvKysqCrtVfaGioq6l1jcfserQ4hIiIXkS2i0he3e8rReR7EdlRdxvYxHbOBm0+9rpiIiJqF5mZOJeXh9J583A4Lw+lmZnw+mq6Xbp0cd1+++2V06dPj7/vvvuaPCD16NGjxo0bN3YEgHfeeSdi6NChVQMHDrx4+PBhY2FhYTAA/PWvf+12xx13tDgYhYWFOX/88UcdADidThw4cMA4atSoyuXLl/9QWVmpr6iouOJbHhUVFfpevXrV6PV6LF++vFt9mCgtLTXGxMTUzJkz5+SkSZNObtu2LfTo0aMGp9OJqVOnnl20aNHh3bt3hwJAXFzcpS1btnQEgNWrV4e38uVqkfT09KpVq1ZFAMAHH3zQ+ccff7zq2yrr16/vbDabzx87dmzX4cOHdx85cmT3yJEjz7zzzjtdFyxYUF5cXFxUXFxcFB8fX9OS/hobPnx41SeffNK1qqpKzpw5o9uwYUPX63kurZkJeRjA3kb3zVNKDay77WhiuwsN2oy+niKJiOjGyMzEuZwcHGuLAFJvwoQJp0tKSkKa2hUDALfeeuvFN954o5vJZLKcOXPGMHfu3PLQ0FD16quvlo0fP76vyWSy6HQ6zJ07t7yljztx4sTTubm5PRMTEy2FhYXB2dnZfUwmkyU5OdkyY8aM442/5fHII4+cWLNmTbfU1FRzaWlph5CQEBcA/P3vf+9ksViSEhMTLevXrw+fP3/+8bKysqBhw4YlmM1my7Rp0/o888wzPwDAE088cfwvf/lL1KBBg8wnT55sl+uzPffcc0c+//zzzhaLJfGTTz7pEhUVVdO1a9crnss777wTMXr06Ct2kYwdO/bM2rVrrzr2pCX9NTZs2LDzY8aMOZ2cnJyUlZXVNz09vepa7ZsiLZnWEZFYAG8BWAzgMaVUloisBJCnlHq/mW2rlFJhLS0oLS1NFRQUtLQ5EREBEJGtSqm0hvft3LmzLDU19aRWNbVUSUmJMSsrq/++ffv2aF2LL7hw4YIYDAYVFBSEjRs3dpw9e/YtxcXFRTdLf43t3LkzMjU1Nd7TupamtKUA5gNXHaS0WESeApAP4Aml1CUP23YQkQIAtQCeU0p91LiBiDwI4EEA6N3b63PYEBER+a39+/cb77///r4ulwtBQUHqtddeK7uZ+muNZmdCRCQLwD1KqVkiYgMwt24mpBeAYwCMAFYAOKCUesbD9tFKqSMiciuAzwFkKKUONPV4nAkhImo9X54JIf92rZmQlhwT8lMAo0WkDMC7AIaLyNtKqaPK7RKANwGke9pYKXWkbvkvAHYAg1r9DIiIiMjvNBtClFILlFKxSql4ABMAfK6UmlQ3EwJxf0fnXgCFjbcVkXARCa77ORLuQNNm+5mIiIjId3lz5O5qEYkCIAB2AJgJACKSBmCmUmo6gEQAr4mIC+7A85xSiiGEiIiIWhdClFJ2uHepQCk1vIk2BQCm1/28GUCKVxUSERGRX+IZU4mIqN38+9//NmRlZd0aFxeX3Ldv36S77rqr365du67rYmf1Bg0aZPZ0f8OLv3lj9erVXX73u9/1vFabkpIS46uvvurVlWnJu90xRERETXK5XBg9enS/7OzsU3l5ef8CgM2bN4ccOXIkqP56MK1RW1sLg8GA7du3F7d9tf8xceLECgAV12qzb9++4LVr10bMnDmzyROwUfM4E0JERJdt3Lix44IFC3rWn0LdG3l5eZ0MBoOaP3/+5bOcDh069MLIkSOvOrvmnj17glNTU83JycmJjzzySHRoaOig+j6GDBliGjVqVJ+EhIQkAKhf53K5MGXKlN59+/ZNstls/Zo6Q2l6enrCtGnT4gYNGmTu379/0hdffBEKAMePH9dnZmb2NZlMltTUVPO3334bAgC5ubndpkyZ0htwz65MnTo1btCgQebY2NiU+pmWJ598MqagoCDMbDZb/vjHP3YvKCjokJKSkmg2my0mk8mye/dur2Z7AgVDCBERAXAHkKysLNMLL7wQk5WVZfI2iOzatSskNTXV4wXrGps9e3bcrFmzThQWFu6Njo6+4spyu3bt6rhkyZLDBw4cuOKMqqtWreq6f//+4JKSkj0rV648uG3btibPzn3+/Hnd9u3bi3Nzcw8++OCDfQBg/vz50ampqedLS0uLnn322cMPPPBAH0/bHj9+PKigoKB4/fr1+55++ukYAFi8ePHhtLS0quLi4qKnn376xLJly6JmzZp1vLi4uGjXrl17+/TpU+2pL7oSQwgREQEA8vPzO9XU1OhcLhdqa2t1+fn5jc+S3W62b98eNm3atNMAMH369FMN1w0YMOCc2Wy+6k39yy+/7HT//fefNhgMiI+Pr7FarU1e3C47O/s04L4iblVVle7kyZP6LVu2dPrNb35zCgBGjx5defbsWcOpU6euunjb6NGjz+r1egwePPjiqVOngjz1b7Vaz7300ku9nnzyyZ779u0zhoWFtewSwAGOIYSIiAAAGRkZlUFBQS69Xg+DweDKyMho8RVrPUlJSbmwc+fOUE/rfvvb38bUX06+uX7a4nLyjduJCDydMVxErrqzQ4cOl+9r6izjM2fOPL1+/fr9ISEhrrvvvtv08ccf37AA58sYQnydwwHk5LiXREReyMzMPJeXl1c6b968w3l5eaWZmZleXUl31KhRldXV1fLSSy9F1t/35Zdfhn7yySdhy5YtO1x/OXkAGDhwYNXKlSvDAeCNN95o0bdO7rrrrsr33nsvora2FgcPHgz65z//2eQb/5o1a8IB4O9//3tYp06dnN26dXPefvvtlW+++WY3wH3sSXh4eG1ERESTgaehLl26OKuqqi7Pmt+s9mwAAAy/SURBVBQVFRkTExMv/f73vz8xYsSIszt27AhpST+Bjt+O8WUOB5CRAVRXA0YjkJ8PWK1aV0VtyeEA7HbAZuPY+qubbIwzMzPPeRs+6ul0Onz88ccHZs2aFbd06dKewcHBKjY29tKyZcsONW67bNmyQxMnTuyTm5vbc8SIEWfDwsKueSl5AJg8efLZ/Pz8zgkJCUl9+vS5mJ6e3uTMTXh4uHPQoEHmqqoq/YoVK74HgOeff/5IdnZ2vMlksoSEhLhWrlz5fUufW3p6+gWDwaASEhIs2dnZJy9evKh77733uhkMBhUVFVWTk5NzpKV9BbJmL2B3o/ECdq2QkwMsXAg4nYBeDzz7LLBggdZVUVthyPR/bTjGvn4Bu8rKSl3Hjh1dOp0OK1asCF+7dm1Efn5+kxc7bY309PSEF1988dCdd97ZooNkqW1d6wJ2nAnxZTab+w9X/R8wm03riqgt2e3usXU63Uu7nSHE33CML/vmm29CH3744d5KKXTu3Nm5cuXKMq1rovbHEOLLrFb3J6ebaCqX2hBDpv/jGF82cuTIqpKSkna5ttiWLVtK2qNf8h5DiK+zWhk+/BVDpv9r/zF2uVwu0el0N9d+dwoYLpdLADR5sC9DCNHNjCHT/7XvGBeWl5dboqKiKhhE6EZzuVxSXl7eBUBhU20YQoiI/FRtbe30Y8eOvX7s2LFk8JQMdOO5ABTW1tZOb6oBQwgRkZ8aPHjwCQCjta6DqClMxkRERKQJhhAiIiLSBEMIERERaYIhhIiIiDTBEEJERESaYAghIiIiTTCEEBERkSYYQoiIiEgTDCFERESkCYYQIiIi0gRDCBEREWmCIYSIiIg0wRBCREREmmAIISIiIk0whBAREZEmGEKIiIhIEwwhREREpAmGECIiItJEi0OIiOhFZLuI5NX9vlJEvheRHXW3gdfYtrOIHBaRV9qiaCIiIvJ9hla0fRjAXgCdG9w3Tyn1fgu2fRbAl60pjIiIiPxbi2ZCRCQWwC8AvN7aBxCRwQB6APhHa7clIiIi/9XS3TFLAcwH4Gp0/2IR2SUiL4tIcOONREQH4CUA87wrk4iIiPxNsyFERLIAnFBKbW20agEAM4CfAIgA8LiHzWcB+D+l1KFmHuNBESkQkYLy8vKWVU5EREQ+rSXHhPwUwGgRuQdABwCdReRtpdSkuvWXRORNAHM9bGsFcIeIzAIQBsAoIlVKqScaNlJKrQCwAgDS0tLUdT4XIiIi8iHNhhCl1AK4Zz0gIjYAc5VSk0Skl1LqqIgIgHsBFHrYdmL9zyIyFUBa4wBCREREgcmb84SsFpHdAHYDiASwCABEJE1EWn0AKxEREQUWUerm2vuRlpamCgoKtC6DiMiniMhWpVSa1nUQtQbPmEpERESaYAghIiIiTTCEEBERkSb8KoQ4HEBOjnsZCBwOB3JycuAIlCdMRER+pTXXjrmpORxARgZQXQ0YjUB+PmC1al1V+3E4HMjIyEB1dTWMRiPy8/Nh9ecnXMfhAOx2wGbz7/Gt53A4YLfbYbPZOL5E5Hf8JoTY7e4A4nS6l3a7f/8Rs9vtqK6uhtPpRHV1Nex2u9+/STFo+nfQDLTxrRdoQZOoIb/ZHWOzuf9w6fXupc2mdUXty2azwWg0Qq/Xw2g0wubvTxieg6Y/8xQ0/VmgjS/wn6C5cOFCZGRkcNcqBRy/mQmxWt2fnAJlKtdqtSI/Pz+gPkHVB836T8r+nrvqg2b9TIi/B81AG18gMGc0iRriycrIpwTaMQOBNlUfiOPbVrvceLIy8kUMIUREGmqroMkQQr7Ib3bHEBH5IqvVGhCzXESe+M2BqURERORbGEKIiIhIEwwhREREpAmGECIiItIEQwgRERFpgiGEiIiINMEQQkRERJpgCCEiIiJNMIQQERGRJhhCiIiISBMMIURERKQJhhAiIiLSBEMIERERaYIhhIiIiDTBEEJERESaYAghIiIiTTCEEBERkSYYQoiIiEgTDCFERESkCYYQIiIi0gRDCBEREWmCIYSIiIg0wRBCREREmmhxCBERvYhsF5G8ut9Xisj3IrKj7jbQwza3iMjWuvV7RGRmWxZPREREvsvQirYPA9gLoHOD++Yppd6/xjZHAQxVSl0SkTAAhSLysVLqyHXUSkRERH6kRTMhIhIL4BcAXm9N50qpaqXUpbpfg1v6eEREROT/WhoKlgKYD8DV6P7FIrJLRF4WkWBPG4pInIjsAnAIwPOeZkFE5EERKRCRgvLy8tbUT0RERD6q2RAiIlkATiiltjZatQCAGcBPAEQAeNzT9kqpQ0qpAQD6AXhARHp4aLNCKZWmlEqLiopq7XMgIiIiH9SSmZCfAhgtImUA3gUwXETeVkodVW6XALwJIP1andTNgOwBcIeXNRMREZEfaDaEKKUWKKVilVLxACYA+FwpNUlEegGAiAiAewEUNt5WRGJFJKTu53C4A01JG9ZPREREPqo1345pbLWIRAEQADsAzAQAEUkDMFMpNR1AIoCXRETVtXtRKbXby5qJiIjID4hSSusarpCWlqYKCgq0LoOIyKeIyFalVJrWdRC1Br8y68scDiAnx70k/8QxJiI/5s3uGNKSwwFkZADV1YDRCOTnA1ar1lVRW+IY+z+HA7DbAZuNY0sBiTMhvspud785OZ3upd2udUXU1jjG/q0+ZC5c6F5ytosCEEOIr7LZ3J+O9Xr30mbTuiJqaxxj/8aQScTdMT7LanVPz3Mq139xjP1bfcis393GkEkBiN+OISLSShseE8Jvx5Av4kwIEZFWrFbOcFFA4zEhREREpAmGECIiItIEQwgRERFpgiGEiIiINMEQQkRERJpgCCEiIiJN3HTnCRGRcgAHr9EkEsDJG1ROe/Dl+n25dsC362ft2vGV+m9RSkVpXQRRa9x0IaQ5IlLgyyfk8eX6fbl2wLfrZ+3a8fX6iW5m3B1DREREmmAIISIiIk34YghZoXUBXvLl+n25dsC362ft2vH1+oluWj53TAgRERH5B1+cCSEiIiI/wBBCREREmrgpQ4iILBGRYhHZJSIfikjXRut7i0iViMxtYvtNIrKj7nZERD66MZVffnxv6xcRWSwipSKyV0T+98ZU3ia1rxSR7xu8/gNvTOXe196g3TIRqWrfaj0+rrev/V9EZGfd9u+LSNiNqbxNal8tIiUiUigib4hI0I2p/PLje1v/bBHZLyJKRCJvTNVEvu+mDCEANgBIVkoNAFAKYEGj9S8D+LSpjZVSdyilBiqlBgJwAPig3Sr1zKv6AUwFEAfArJRKBPBuexTZBG9rB4B59a+/UmpHexTZBK9rF5E0AF2v1aYdeVv/o0qp1Lrt/w1gdvuU6ZG3ta8GYAaQAiAEwPT2KPIavK3/GwCZuPaJFomokZsyhCil/qGUqq379Z8AYuvXici9AP4FYE9z/YhIJwDDAdzQmZA2qP//AXhGKeWq6+9Ee9XaWFu99lrwtnYR0QNYAmB+e9bZFG/rV0r9WNdW4H4jv2FHnbdB7f+n6gDY0nD7G6EN6t+ulCpr1yKJ/NBNGUIamYa6TyAi0hHA4wD+2MJtxwDIr//jrJHrqb8vgF+KSIGIfCoi/du5xqZc72u/uG5a+2URCW7PAq/hemqfDeBjpdTRdq6tJa7rtReRNwEcg3tWYVl7FngN1/1/tm43zGQAn7Vbdc3z5m8OEbWCQasHFpGNAHp6WPWkUmp9XZsnAdTCPVULuP8QvKyUqnJ/2GvWrwC83gblXqWd6w8GcFEplSYi9wF4A8AdPlL7ArjfBI1wn1/hcQDPtFHp7Va7iEQDGA/A1la1NvE47frvXin167oZnWUAfgngTV+pvc5yAF8ppTa1QclXuEH1E1FrKKVuyhuAB+A+niO0wX2bAJTV3c4COA1gdhPbdwNwCkAHX6sfQDGA+LqfBUCFr9TeqB8bgDxfqB3AL+AOT/XtXAD2+9K/m0b93OUrr32Dtk/DvetUd6Nf97Z67evaRWpRP2+8+eJN8wI8FgWMBFAEIOoabf4AYO411s8E8JYv1g/gOQDT6n62AfjOh2rvVbcUAEsBPOcrtTdqV+VL/27qXu9+DX5+EcCLvlB73brpADYDCLnRr3tb/tthCOGNt9bdbtZjQl4B0AnAhrqveb7a3AYi8n91U+r1JgBY014FNsPb+p8DMFZEdgPIwY39poC3ta+uq3s33JdAX9R+pV6lLf7daMmb+gXAWw1e+15ow91gLeDta/8qgB4AHHXbP9WOtXriVf0i8r8i8gPcB7TuEpF22Q1M5G942nYiIiLSxM06E0JERER+jiGEiIiINMEQQkRERJpgCCEiIiJNMIQQERGRJhhCiIiISBMMIURERKSJ/w8sm8ViarZoKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "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(fieldsetA.U.grid.lon[:nind, :nind], fieldsetA.U.grid.lat[:nind, :nind], '.r', label='U points assuming A-grid')\n",
    "ax2 = ax.plot(fieldsetA.V.grid.lon[:nind, :nind], fieldsetA.V.grid.lat[:nind, :nind], '.b', label='V points assuming A-grid')\n",
    "\n",
    "ax3 = ax.plot(fieldsetC.U.grid.lon[:nind, :nind], fieldsetC.U.grid.lat[:nind, :nind], '.k', label='C-grid points')\n",
    "\n",
    "ax.legend(handles=[ax1[0], ax2[0], ax3[0]], loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further information about the NEMO C-grids is available in [the NEMO 3D tutorial](https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_3D.ipynb)."
   ]
  },
  {
   "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`."
   ]
  },
  {
   "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."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}