{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NestedFields\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some applications, you may have access to different fields that each cover only part of the region of interest. Then, you would like to combine them all together. You may also have a field covering the entire region and another one only covering part of it, but with a higher resolution. The set of those fields form what we call nested fields.\n",
    "\n",
    "It is possible to combine all those fields with kernels, either with different if/else statements depending on particle position.\n",
    "\n",
    "However, an easier way to work with nested fields in Parcels is to combine all those fields into one `NestedField` object. The Parcels code will then try to successively interpolate the different fields.\n",
    "\n",
    "For each Particle, the algorithm is the following:\n",
    "\n",
    "1. Interpolate the particle onto the first `Field` in the `NestedFields` list.\n",
    "\n",
    "2. If the interpolation succeeds or if an error other than `ErrorOutOfBounds` is thrown, the function is stopped.\n",
    "\n",
    "3. If an `ErrorOutOfBounds` is thrown, try step 1) again with the next `Field` in the `NestedFields` list\n",
    "\n",
    "4. If interpolation on the last `Field` in the `NestedFields` list also returns an `ErrorOutOfBounds`, then the Particle is flagged as OutOfBounds.\n",
    "\n",
    "This algorithm means that **the order of the fields in the 'NestedField' matters**. In particular, the smallest/finest resolution fields have to be listed _before_ the larger/coarser resolution fields.\n",
    "\n",
    "Note also that `pset.execute()` can be _much_ slower on `NestedField` objects than on normal `Fields`. This is because the handling of the `ErrorOutOfBounds` (step 3) happens outside the fast inner-JIT-loop in C; but instead is delt with in the slower python loop around it. This means that especially in cases where particles often move from one nest to another, simulations can become very slow.\n",
    "\n",
    "This tutorial shows how to use these `NestedField` with a very idealised example.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "from parcels import AdvectionRK4, Field, FieldSet, JITParticle, NestedField, ParticleSet"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First define a zonal and meridional velocity field defined on a high resolution (dx = 100m) 2kmx2km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 _ cos(lon / 200 _ pi / 2) m/s.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 21\n",
    "lon = np.linspace(0.0, 2e3, dim, dtype=np.float32)\n",
    "lat = np.linspace(0.0, 2e3, dim, dtype=np.float32)\n",
    "lon_g, lat_g = np.meshgrid(lon, lat)\n",
    "V1_data = np.cos(lon_g / 200 * np.pi / 2)\n",
    "U1 = Field(\"U1\", np.ones((dim, dim), dtype=np.float32), lon=lon, lat=lat)\n",
    "V1 = Field(\"V1\", V1_data, grid=U1.grid)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define the same velocity field on a low resolution (dx = 2km) 20kmx4km grid.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "xdim = 11\n",
    "ydim = 3\n",
    "lon = np.linspace(-2e3, 18e3, xdim, dtype=np.float32)\n",
    "lat = np.linspace(-1e3, 3e3, ydim, dtype=np.float32)\n",
    "lon_g, lat_g = np.meshgrid(lon, lat)\n",
    "V2_data = np.cos(lon_g / 200 * np.pi / 2)\n",
    "U2 = Field(\"U2\", np.ones((ydim, xdim), dtype=np.float32), lon=lon, lat=lat)\n",
    "V2 = Field(\"V2\", V2_data, grid=U2.grid)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now combine those fields into a `NestedField` and create the fieldset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = NestedField(\"U\", [U1, U2])\n",
    "V = NestedField(\"V\", [V1, V2])\n",
    "fieldset = FieldSet(U, V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Output files are stored in NestedFieldParticle.zarr.\n",
      "100%|██████████| 14000.0/14000.0 [00:01<00:00, 9625.92it/s] \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAGdCAYAAAAWp6lMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHp0lEQVR4nO3df1xUdb4/8NcoMALBCUQYJlFZ19KCrNBVzMRfgW5ErX1XXbpc3TWy9Vcuekvq3kfUTTEr625uZq5ppa7uvWW/NBTLHxGgSFL4MzMUTEbUhRlRHFA/3z9sTgwzA6jMr8+8no/HPB7OmQ+H83nPkXnN53zOORohhAARERGRxDq5ewOIiIiInI2Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikp6fuzfAU1y5cgUnT55ESEgINBqNuzeHiIiI2kEIgXPnzkGv16NTJ8fjOAw8Pzt58iRiYmLcvRlERER0HaqqqtC9e3eHrzPw/CwkJATA1YKFhoa6eWuIiIioPUwmE2JiYtTPcUcYeH5mOYwVGhrKwENERORl2pqOwknLREREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9pwaepUuX4s4771SvXpyYmIjPP/9cfV0IgZycHOj1egQGBmL48OHYv3+/1TrMZjNmzpyJiIgIBAcHIy0tDSdOnLBqU1tbi4yMDCiKAkVRkJGRgbq6Omd2jYiIiLyIUwNP9+7dsXDhQuzZswd79uzByJEj8dBDD6mhZtGiRVi8eDGWLFmCkpIS6HQ63H///Th37py6jtmzZ2PDhg1Yt24dCgoKUF9fj9TUVFy+fFltk56ejrKyMuTl5SEvLw9lZWXIyMhwZteIiIjIi2iEEMKVvzA8PBwvv/wy/vSnP0Gv12P27Nl4+umnAVwdzYmKisJLL72EqVOnwmg0olu3bnj//fcxYcIEAL/c1XzTpk1ISUnBwYMHcfvtt6O4uBiDBg0CABQXFyMxMRGHDh3Cbbfd1q7tMplMUBQFRqOxQ++lJYTAhaYLHbY+IiIibxTkH9Tm/a6uR3s/v11289DLly/jf//3f3H+/HkkJiaioqICBoMBycnJahutVoukpCQUFhZi6tSpKC0tRVNTk1UbvV6PuLg4FBYWIiUlBUVFRVAURQ07ADB48GAoioLCwkKHgcdsNsNsNqvPTSZTh/dZCIGhK4eisKqww9dNRETkTe6NuRdf/fErp4Se9nD6pOXy8nLcdNNN0Gq1eOKJJ7BhwwbcfvvtMBgMAICoqCir9lFRUeprBoMBAQEBCAsLa7VNZGSkze+NjIxU29iTm5urzvlRFAUxMTE31E97LjRdYNghIiIC8HXV12494uH0EZ7bbrsNZWVlqKurwwcffIBJkyZhx44d6ustk54Qos3017KNvfZtrSc7OxtZWVnqc5PJ5JTQY3Fq7ikE+wc7bf1ERESe6HzTeUS9EtV2QydzeuAJCAjAr3/9awDAgAEDUFJSgv/5n/9R5+0YDAZER0er7WtqatRRH51Oh8bGRtTW1lqN8tTU1GDIkCFqm1OnTtn83tOnT9uMHjWn1Wqh1WpvvIPtFOwfjOAABh4iIiJ3cPl1eIQQMJvNiI2NhU6nQ35+vvpaY2MjduzYoYaZhIQE+Pv7W7Wprq7Gvn371DaJiYkwGo3YvXu32mbXrl0wGo1qGyIiIvJtTh3heeaZZzB27FjExMTg3LlzWLduHbZv3468vDxoNBrMnj0bCxYsQJ8+fdCnTx8sWLAAQUFBSE9PBwAoioIpU6Zgzpw56Nq1K8LDwzF37lzEx8dj9OjRAIB+/fphzJgxyMzMxLJlywAAjz/+OFJTU9t9hhYRERHJzamB59SpU8jIyEB1dTUURcGdd96JvLw83H///QCAp556Cg0NDZg2bRpqa2sxaNAgbNmyBSEhIeo6XnvtNfj5+WH8+PFoaGjAqFGjsGrVKnTu3Flts2bNGsyaNUs9mystLQ1LlixxZteIiIjIi7j8OjyeyhnX4TnfeB435d4EAKjPruccHiIi8jnO/ixs7+c376VFRERE0mPgISIiIukx8BAREZH0GHiIiIhIegw8REREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUmPgYeIiIikx8BDRERE0mPgISIiIukx8BAREZH0GHiIiIhIegw8REREJD0GHiIiIpIeAw8RERFJz6mBJzc3FwMHDkRISAgiIyPx8MMP4/Dhw1ZtJk+eDI1GY/UYPHiwVRuz2YyZM2ciIiICwcHBSEtLw4kTJ6za1NbWIiMjA4qiQFEUZGRkoK6uzpndIyIiIi/h1MCzY8cOTJ8+HcXFxcjPz8elS5eQnJyM8+fPW7UbM2YMqqur1cemTZusXp89ezY2bNiAdevWoaCgAPX19UhNTcXly5fVNunp6SgrK0NeXh7y8vJQVlaGjIwMZ3aPiIiIvISfM1eel5dn9XzlypWIjIxEaWkphg0bpi7XarXQ6XR212E0GrFixQq8//77GD16NABg9erViImJwdatW5GSkoKDBw8iLy8PxcXFGDRoEABg+fLlSExMxOHDh3Hbbbc5qYdERETkDVw6h8doNAIAwsPDrZZv374dkZGRuPXWW5GZmYmamhr1tdLSUjQ1NSE5OVldptfrERcXh8LCQgBAUVERFEVRww4ADB48GIqiqG1aMpvNMJlMVg8iIiKSk8sCjxACWVlZGDp0KOLi4tTlY8eOxZo1a/Dll1/i1VdfRUlJCUaOHAmz2QwAMBgMCAgIQFhYmNX6oqKiYDAY1DaRkZE2vzMyMlJt01Jubq4630dRFMTExHRUV4mIiMjDOPWQVnMzZszAd999h4KCAqvlEyZMUP8dFxeHAQMGoGfPnti4cSPGjRvncH1CCGg0GvV58387atNcdnY2srKy1Ocmk4mhh4iISFIuGeGZOXMmPvnkE2zbtg3du3dvtW10dDR69uyJI0eOAAB0Oh0aGxtRW1tr1a6mpgZRUVFqm1OnTtms6/Tp02qblrRaLUJDQ60eREREJCenBh4hBGbMmIEPP/wQX375JWJjY9v8mbNnz6KqqgrR0dEAgISEBPj7+yM/P19tU11djX379mHIkCEAgMTERBiNRuzevVtts2vXLhiNRrUNERER+S6nHtKaPn061q5di48//hghISHqfBpFURAYGIj6+nrk5OTgkUceQXR0NI4dO4ZnnnkGERER+N3vfqe2nTJlCubMmYOuXbsiPDwcc+fORXx8vHrWVr9+/TBmzBhkZmZi2bJlAIDHH38cqampPEOLiIiInBt4li5dCgAYPny41fKVK1di8uTJ6Ny5M8rLy/Hee++hrq4O0dHRGDFiBNavX4+QkBC1/WuvvQY/Pz+MHz8eDQ0NGDVqFFatWoXOnTurbdasWYNZs2apZ3OlpaVhyZIlzuweEREReQmNEEK4eyM8gclkgqIoMBqNHTaf53zjedyUexMAoD67HsEBwR2yXiIiIm/h7M/C9n5+815aREREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUmPgYeIiIikx8BDRERE0mPgISIiIukx8BAREZH0GHiIiIhIegw8REREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9JwaeHJzczFw4ECEhIQgMjISDz/8MA4fPmzVRgiBnJwc6PV6BAYGYvjw4di/f79VG7PZjJkzZyIiIgLBwcFIS0vDiRMnrNrU1tYiIyMDiqJAURRkZGSgrq7Omd0jIiIiL+HUwLNjxw5Mnz4dxcXFyM/Px6VLl5CcnIzz58+rbRYtWoTFixdjyZIlKCkpgU6nw/33349z586pbWbPno0NGzZg3bp1KCgoQH19PVJTU3H58mW1TXp6OsrKypCXl4e8vDyUlZUhIyPDmd0jIiIiL6ERQghX/bLTp08jMjISO3bswLBhwyCEgF6vx+zZs/H0008DuDqaExUVhZdeeglTp06F0WhEt27d8P7772PChAkAgJMnTyImJgabNm1CSkoKDh48iNtvvx3FxcUYNGgQAKC4uBiJiYk4dOgQbrvttja3zWQyQVEUGI1GhIaGdkh/zzeex025NwEA6rPrERwQ3CHrdTohgMsX3L0V16dzEKDRuHsriIjoZ87+LGzv57dfh/7WNhiNRgBAeHg4AKCiogIGgwHJyclqG61Wi6SkJBQWFmLq1KkoLS1FU1OTVRu9Xo+4uDgUFhYiJSUFRUVFUBRFDTsAMHjwYCiKgsLCQruBx2w2w2w2q89NJlOH99crCQHkDwXOFLp7S65Pt3uB0V8x9BARkRWXTVoWQiArKwtDhw5FXFwcAMBgMAAAoqKirNpGRUWprxkMBgQEBCAsLKzVNpGRkTa/MzIyUm3TUm5urjrfR1EUxMTE3FgHZXH5gveGHQA4/bX3jk4REZHTuGyEZ8aMGfjuu+9QUFBg85qmxbdxIYTNspZatrHXvrX1ZGdnIysrS31uMpkYeloadwrw85LDcJfOAx9Gtd2OiIh8kksCz8yZM/HJJ59g586d6N69u7pcp9MBuDpCEx0drS6vqalRR310Oh0aGxtRW1trNcpTU1ODIUOGqG1OnTpl83tPnz5tM3pkodVqodVqb7xzMvML9p7AQ0RE1AqnHtISQmDGjBn48MMP8eWXXyI2Ntbq9djYWOh0OuTn56vLGhsbsWPHDjXMJCQkwN/f36pNdXU19u3bp7ZJTEyE0WjE7t271Ta7du2C0WhU2xAREZHvcuoIz/Tp07F27Vp8/PHHCAkJUefTKIqCwMBAaDQazJ49GwsWLECfPn3Qp08fLFiwAEFBQUhPT1fbTpkyBXPmzEHXrl0RHh6OuXPnIj4+HqNHjwYA9OvXD2PGjEFmZiaWLVsGAHj88ceRmprarjO0iIiISG5ODTxLly4FAAwfPtxq+cqVKzF58mQAwFNPPYWGhgZMmzYNtbW1GDRoELZs2YKQkBC1/WuvvQY/Pz+MHz8eDQ0NGDVqFFatWoXOnTurbdasWYNZs2apZ3OlpaVhyZIlzuweEREReQmXXofHk/E6PD+7dB7459Vtxvh675nD463bTUQkOU+5Dg/vpUVERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUmPgYeIiIikx8BDRERE0mPgISIiIukx8BAREZH0GHiIiIhIegw8REREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUnPqYFn586dePDBB6HX66HRaPDRRx9ZvT558mRoNBqrx+DBg63amM1mzJw5ExEREQgODkZaWhpOnDhh1aa2thYZGRlQFAWKoiAjIwN1dXXO7BoRERF5EacGnvPnz6N///5YsmSJwzZjxoxBdXW1+ti0aZPV67Nnz8aGDRuwbt06FBQUoL6+Hqmpqbh8+bLaJj09HWVlZcjLy0NeXh7KysqQkZHhtH4RebpqYwMKj55BtbHB3ZsiNdbZdVhrulF+zlz52LFjMXbs2FbbaLVa6HQ6u68ZjUasWLEC77//PkaPHg0AWL16NWJiYrB161akpKTg4MGDyMvLQ3FxMQYNGgQAWL58ORITE3H48GHcdtttHdspIg9VbWzAnmP/whcHa/Bx2UmIn5en/yYGEwbG4HzjZcRGBCNaCXTrdno7e3XWAPjDb2LQTx+KsKAAJPQMY51vkKXOdQ1NOFh9Dut2V+LKzzv1w3dFY/TtOtaZrolTA097bN++HZGRkbj55puRlJSE+fPnIzIyEgBQWlqKpqYmJCcnq+31ej3i4uJQWFiIlJQUFBUVQVEUNewAwODBg6EoCgoLCx0GHrPZDLPZrD43mUxO6iGR8y3beRS5mw7ZfW3t7iqs3V0F4OoH87yxfTE1qbcLt04ejuosALXGFtms83VrbX8GgI/KqvFRWTUA1pnaz62TlseOHYs1a9bgyy+/xKuvvoqSkhKMHDlSDSIGgwEBAQEICwuz+rmoqCgYDAa1jSUgNRcZGam2sSc3N1ed86MoCmJiYjqwZ0Su8/LmQ61+ODQnAOR+fggvb25fe/rFtdQZYJ2vF+tMzuLWEZ4JEyao/46Li8OAAQPQs2dPbNy4EePGjXP4c0IIaDQa9Xnzfztq01J2djaysrLU5yaTiaGHvEq1sQF//eII/tFiZKE9/rbtKADgP1L6dvRmSenlzYfUml0L1vnasM7kTG4/pNVcdHQ0evbsiSNHjgAAdDodGhsbUVtbazXKU1NTgyFDhqhtTp06ZbOu06dPIyoqyuHv0mq10Gq1HdwDItdYX1KJeR+Uq/N0rsffth1FaKA/pg7j4YDWLNtx9Lo+hC1Y5/ZhncnZPOo6PGfPnkVVVRWio6MBAAkJCfD390d+fr7aprq6Gvv27VMDT2JiIoxGI3bv3q222bVrF4xGo9qGSCbVxgZkf3hjYcdi4aZDPOulFdXGBuR+fuOHS1jn1lUbG7CQdSYnc+oIT319PX744Qf1eUVFBcrKyhAeHo7w8HDk5OTgkUceQXR0NI4dO4ZnnnkGERER+N3vfgcAUBQFU6ZMwZw5c9C1a1eEh4dj7ty5iI+PV8/a6tevH8aMGYPMzEwsW7YMAPD4448jNTWVZ2iRlN748oh6too9KbdHIe0uPbqHBaLqXw3QaICSin9hVdFxm7YCwBtf/IAF4+Kdt8FeqtrYgOc+2u/wdUud7+l5dfS59FgtPv32JDYfsB1xZp1b9z9bjzgM8Nda5/mfHcSzqf149hbZcGrg2bNnD0aMGKE+t8yZmTRpEpYuXYry8nK89957qKurQ3R0NEaMGIH169cjJCRE/ZnXXnsNfn5+GD9+PBoaGjBq1CisWrUKnTt3VtusWbMGs2bNUs/mSktLa/XaP0TeatmOo1i7y/6cHXtnYPWPufoh8cCdegR38bN7yGDt7kr0jAjioYBm2jpkOH1Eb5v5Iqn9A5HaX+9wHgrrbN+yHUexrsT+Pn09df6svBqb9lUjd1w8Jgzs4ZRtJu+kEUJ0xMi41zOZTFAUBUajEaGhoR2yzvON53FT7k0AgPrsegQHBHfIep3q0nngn1e3GePrAT8v2GbAe7f7GlQbGzAk90u7H8Kp8dHt+lb7zIbv7AYmDYDC7JH8Voyrdb534ZcOR9HSf9OjzZGa1ur80fQhahD1dd9W1eKhvxXafe1G6gxcna/xNfdpj+Dsz8L2fn571BweInJsZUGF3bCjAdo9hD9zZB/YO3dRAFhZcOzGNlASKwsqHIadTgBmjvp1m+torc4Pv1mI9SWVN7KJUlhfUomHHYQdDdpf504OTsa9Au7TZI2Bh8gLVBsbsPyrCruvzftt33Z/i41WAjFvrP1Td5d/9aPPT/hsrc6dNEDuI/HtqnVrdRYCyP6g3Kdr3dbE+/bu09FKIHLHxTv8IOM+Tc0x8BB5AUejO+m/6XHNc0KmJvVG+iDba05xlMdxnVPjo/H1vJHXNCdkalJvZP+2r92RHl8ffXA0iqbBz1dOvoZ9esLAHvg6eyRS46NtXuM+Tc0x8BB5OEejDu0d9rfH0SGX5V/9iG+raq9rnd6utTpf71k/U4f1xkfT7V8ew1dHH1qr80fTh1zXbSKilUA8m9rP4T7ti3UmWww8RB7O0ahD5n2/uu4JmdFKIDLvi7VZ7stzTJxRZ+DqmXKPO6i1L44+tFbnG5nM3do+7Yt1JlsMPEQerLVvw38c2uuG1v3HobF2J3wKATzz4T6f+lbszDoDV2ttb/Th7wW+NfrAOpM7MfAQebCKM+edMuoAtD7h87IQOHbmwg2t35s4s86A49GHKwKsM1hncg0GHiIPtreyzmZZJ3TMt2Hg6oTPDdOHwN59dr/7yfZ3y6q8ymizrCPrDDgeUfOpOp9gncl9GHiIPFS1sQGvbD5ss/zpazgNvT36x4Rh3hjbU6hf8pH7ElUbG7Awz/Y+Th1d52glEE+zzjbLWWdyFQYeIg/1joPJnXfecnOH/6747orNMl85ddrRJFrWuWOtLKiAvev6s87kKgw8RB7I0eTOThqgV0RQh/++2Ihgn5zsyTq7ButMnoCBh8gDlR63fy2cx4Z2zOTOlnx1sqejSbSsc8dinckTMPAQeZj1JZWYuXavzfKOntzZki9O9vzOBZOVW/LFOrtisnJLvlhnah0DD5EHcXSPoWu5j9P1cjTZc9Hnh6U8DFBtbMBLm50/ibYln6yzCyYrt+Rrdaa2MfAQeZCKM+ft3mPorxPvvqb7OF0ve5M9Zb0mT+nxWpdNom3Jl+rsaJ9mncnVGHiIPIi9of/OGg0Sel3/JfevRWxEsE8cBnB02LCzRuOUSbQt+UqdAfvXOGKdyR0YeIg8hKOh/6fG3ubUof/mfOEwQGuHDReMi3NJrX2hzsDP196xc9jQVfu0r9SZ2oeBh8hDuHPovznZDwO4+7Chhex1Bn4+O8vN+7SjOpces38mJMmLgYfIQzg6nOWKof/mZD8M4O7Dhhay1xmwfxacq/dpR3WetW4v1pdUumw7yP0YeIg8gCcczrKQ+TAA6+w6js6Cc3WtHd0k94oAnvlwnxS1pvZh4CHyAJ5yOMtC1sMtrLPreMLhLIsJA3vgr+l32yyXpdbUPgw8RB7A3qXw3XE4q/n2tDwM0AnOuQ2AK8VGBNvcGd4T6xwU4P1/mrv42fbBnbVO6Bkm5T5N7ef9/6vIJaqNDSg8euaGh389bT2eYuf3p63OGtK48IwheyyHAZp/Pghc3U5vtvP701ajDp5S5+YfxFcA/O7NQq+eX7K+pBKPvFVktayzRuMRtW5Ohn2a2s/P3RtAnm/ZzqNY+PkhCAFoAMwb2xdTk3p71nru1V3zejyF5TTp5jQCGHZrNzdt0VXDbu0GjQZqQBC4Oudh2K3d3PahdSM8tc4TBvZAX10IHvpbobrMMr/EG2utnvbfLFh2AvDhtET0j3HtxPCWht3aDRpA/XLh7fs0XRuO8FCrlu04itxNh6w+9HI/P4RlO4961HpWFPx4TevxJPbmlVyB+29yaG+7vPl0Xk+tMwCcb7xss8xb55c4qvOFxitu2Z7m7N3E1FvrTNeOgYccMhgbsPBz27MsAGDhpkPtPpxU7YL1vLr5+3atwxN5yunoLcl2Oq+n1hmQ6xR11pk8FQMPOfR+0TGbb0MWAsDKgmPtWs/KggqXrMcbedJp0i3JdDqvJ9cZkOcUddaZPBkDDzm08uvjrb7+94If2/wjUW1swPKvKlyyHm/kaadJtyTL6byeXmdAjlPUWWfyZAw85FBboyZXRNvzH+wdM3fWeryRp52Obk9CzzCbbfS203m9oc4yXArA0077t0eGOtP1YeCh69aePxKxEcEuW4838rTT0e2JVgKRk3aH1TJvO53XW+rs7adNe9pp//bIUGe6Pgw81C4aDfDIPbdc83VZPvu2ukPWk7fP0OZ6vI2nniZtT/IdUVbPLafzesO8B2+qc8ttYp2dw3J6uoU31Zmun1MDz86dO/Hggw9Cr9dDo9Hgo48+snpdCIGcnBzo9XoEBgZi+PDh2L9/v1Ubs9mMmTNnIiIiAsHBwUhLS8OJEyes2tTW1iIjIwOKokBRFGRkZKCurs6ZXfM5GgH8e2JPq+Hqtv5IVBsbsGDTwQ5ZzwufHmhzPRYGL/mj5cmnSbdUcea8zTJvmffAOruGt9WZp6f7HqcGnvPnz6N///5YsmSJ3dcXLVqExYsXY8mSJSgpKYFOp8P999+Pc+fOqW1mz56NDRs2YN26dSgoKEB9fT1SU1Nx+fIv161IT09HWVkZ8vLykJeXh7KyMmRkZDizaz7nCoCSY7XXdF0We39UnL0eANhbWWe/Ex7G3mE6T5vvYOHN8x68Yf6OhdfX2cPn71jYq7Onbit1HKcGnrFjx+LFF1/EuHHjbF4TQuD111/Hs88+i3HjxiEuLg7vvvsuLly4gLVr1wIAjEYjVqxYgVdffRWjR4/G3XffjdWrV6O8vBxbt24FABw8eBB5eXn4+9//jsTERCQmJmL58uX47LPPcPjwYWd2T0r/V1pld3lnjQYDe9neiwZwfF0We8HjetbzbdW1rWfu/37rFdeJaXkYzxPnO1h487wHb5i/Y+H1dfbw+TsWljp3/jmhufu2F+QabpvDU1FRAYPBgOTkZHWZVqtFUlISCguvXmK9tLQUTU1NVm30ej3i4uLUNkVFRVAUBYMGDVLbDB48GIqiqG3sMZvNMJlMVg9fV21swHMf77dZ3unnP1z9Y8LafV2WamMDXtliGzifGnvbNa9n0eYbX4+n8ab5DhbeOL/EW+vsbfNLvLHOEwb2QMG8EfhH5mAUzBuBCQN7uHuTyMncFngMhquTUKOirCdDRkVFqa8ZDAYEBAQgLCys1TaRkZE264+MjFTb2JObm6vO+VEUBTExMTfUHxk4uobGXyferf4xaO91WSrOnLf6tmdhuR6Hq9fjabxpvoOFN84v8dY6e9v8Em+sM3B1pCexd1eO7PgIt5+lpWlx0FcIYbOspZZt7LVvaz3Z2dkwGo3qo6rK/qEcX+LouHZCL+vAmdAzzOZYfct5BkH+trtWy2Pk7VlPcEDn61qPpx+P96b5OxbeOO/Bm+bvWLDORM7htsCj0129u3XLUZiamhp11Een06GxsRG1tbWttjl16pTN+k+fPm0zetScVqtFaGio1cPXRSuBeP6hX6634ui4drQSiP9+KM5qWfN5ButLKjFuaZHV6/bWFa0EYv7Dra/nd29aH5Z0tJ7m14nRAHhqjGdcyt6RiJu0COj8y0eEN8whsDe/xNPrHK0EQqd0UZ97U52bBwhvqHP/mF+uYOwNdSbf47bAExsbC51Oh/z8fHVZY2MjduzYgSFDhgAAEhIS4O/vb9Wmuroa+/btU9skJibCaDRi9+7daptdu3bBaDSqbaj9/l/CL4f2ts4Z5vC49qh+1ocRLfMMvq2qRfaH5VbD250AfDgt0e66RvTtmPWMHxBjtY6X8g559MTlVzYfRuPlq52zBDRvmEMwYWAPDL/tl3kZnl7ndwoqUG28CMD76vz7Ad3V555e53W7K1FW9ctNQ72lzuRbnBp46uvrUVZWhrKyMgBXJyqXlZWhsrISGo0Gs2fPxoIFC7Bhwwbs27cPkydPRlBQENLT0wEAiqJgypQpmDNnDr744gvs3bsX//Zv/4b4+HiMHj0aANCvXz+MGTMGmZmZKC4uRnFxMTIzM5GamorbbrvNmd2Tnq6Vb2eO5nPYO1X8CoALjVecup6W197x5InL1cYGvL3zR/W5ALAozztuXlhtbMCOZmcMeXqd//uzX67f5G11/r/SX6435ul1fmaD9YRlb6kz+RY/Z658z549GDFihPo8KysLADBp0iSsWrUKTz31FBoaGjBt2jTU1tZi0KBB2LJlC0JCQtSfee211+Dn54fx48ejoaEBo0aNwqpVq9C58y9zO9asWYNZs2apZ3OlpaU5vPYPdQzLPIOWIzADe12dT9N8onFrx/I7aj3Hzp6HrsUyy0RPTxtWLz1e63BSqqdta0v2JpF76ra3NvnX07a1JXuTgD11271pW8m3OXWEZ/jw4RBC2DxWrVoF4Opk45ycHFRXV+PixYvYsWMH4uKs53R06dIFb7zxBs6ePYsLFy7g008/tTmjKjw8HKtXr1ZPL1+9ejVuvvlmZ3bN5zm6Xsh7Rcev6VocHbWeXl3t32vru5/q2uqKS60vqcTMtXttlnvLBE97E2oBz6sz4Pj6Taxzxyo/YbRZ5i11Jt/i9rO0yHvZu17IB9/8ZNWmPdfi6Ij1ODr8tuhzzxlat1yrpOWoQ6c2wpwniVYC8fSYvjbLPanOQOvXb2KdO061sQEv5R2yWe4tdSbfwsBD183eIYOW2nMtjo5ajz2edP2S9lznyBvEd1dslnlSnYG2r9/kDbylzvb2aW+qM/kOBh66bvauvdFSJ03b9wHqqPXY/Tl4zn2I2nudI0/nDfd7as/1mzydN9TZG68nRb6LgYeuW7QSiMz7Yltt89jQX7U5tN1R67HHk+5DtPP701bfhj35XkOt8fT7PbX3+k2ezt71eDypzgCw7VCN1XNv3afJNzDw0A3549BYh6MzGgB/HNrLpeuxtLfwlPsQeeO9hloz7NZuVle39rQ6t/f6TZ7OU+sMXK31sx/ts1rmzfs0yY+Bh25ItBKIeWNtJ1cCwLzf9m33N72OWg8Aj7wPkbfea8iR1k5PdydHdXZ0/SZP19op3+5mbx/w5n2a5MfAQzdsalJvZP+2rzrfoBOA7LF9MXVYb7esxxPnPch2ryFPnV8SGxHsdfdVa42n1hkAeoTbboM315rkx8BDHWLqsN74et5I/CNzML7OHompSdcWUjpyPc3vBwZ4xryHnd+fthp58va5Dp46v2Tn96ev6fpNns5T6wwA7xUet3ru7bUm+THwUIeJVgKR2LvrDf/Bu9H1DP11hNW3YnfPe5Bt/o6Fp80vYZ1dp9rYgLe/+tFqmQy1Jrkx8JB0jp31rHkPss3fsfC0+SW+VufSY7Xu2SAAJRX/slkmQ61Jbgw8JJ1eXT3rsvyyXnrf025/4Gt1nrVur1vuoL6+pBKz1pXZLJeh1iQ3Bh6Sjs6DLssv86X3Pen2B7LXOXdcvM0fa3fcQd3eYUPAu26PQr6LgYek5CmX5Zf90vuO6uzqwy2y13nCwB74a/rdNstdvU/LcnsU8k0MPCQlTzncIuthFgtPOdwie50BIKFnmNv3aUd19rbbo5BvYuAhKXnC4RaZD7NYeMLhFl+oM+D+fdpX6kzyYuAhabn7sJbsh1ks3H24xVfqDLh3n/alOpOcGHhIWu4+rPVdlfyHWSzcebjFFw5nWbhzn/alOpOcGHhIWo4OAby06RC+rXLupFpfG/5vrdbOPNxSbWzAQtbZJXX2pf2Z5MTAQ1KzdwjgCoCH3yx06qTalQUVNjcxBeQe/ndU65UFx5z2O1cWVNjcwBJgnTvayoIKHs4ir8fAQ1JzdAhAOHFSbbWxAcu/qrBZ3knjGTd9dBZ7N0gFgL8X/Mg6dyDWmej6MPCQ1BydRQQ4b7JnxZnzdkd3Hhv6K6mH/6OVQGTeF2uz/Ipwzi0HWGdrrDNR6xh4SHoTBvbAhulDrG7AaOGMyZ7FR8/aLOsE4I9De3X47/I0fxwa67JJtfYmhbPOdR3+u+xNVvaVOpNcGHjIJ/SPCcM8F0z2/LaqFn/98geb5U//tq9PfBt21aRaR5NoWeeOr7O9SeG+UmeSCwMP+QxnT/ZcX1KJh/9WaPc1X5rc6YpJtX//6kefmxTekivq7IuTwkleDDzkMxxN9lz+1Y1P9rTcVNHeh7CvXavEmXUGro6irbDzoe5rk2idXedqYwPe5mRlkggDD/kMR5M9BW78W7Gj03Z98S7Szqxza6NovjaJ1pl1BoB3CmzDDuB7dSZ5MPCQT/nj0NgO/1bs6LRdDYAN04b45F2knVVnR6NovjqJ1hl1Blrfp32xziQHBh7yKa19K37jC9vJxu3xxpdH7H4IZ973K/SP8c27SDurzo5G0XIfiffJUYfW6jz/s4PXFXqqjQ3I+WS/3dcy7+PoDnkvBh7yOY6+Fa/dXYllO49e07qW7TiKtbuqbJbzm7Dr6uyro2gWjur8WXk17l345TVdUXx9SSWG5H6JzftP2bzGfZq8HQMP+RxH34oBYOE1nNZbbWzAws9tT9kF+E0YaLvO7b2fWVt19tVRNIvW6nxFANkflLdrn27tkCHAfZq8HwMP+SRH34qv5ZCLo0NZ/Cb8i9bq/NDfCrFsR9sjPaxz2xxdiBC4eqp6e/ZpR4cMAd+dI0VycXvgycnJgUajsXrodDr1dSEEcnJyoNfrERgYiOHDh2P/fuvjy2azGTNnzkRERASCg4ORlpaGEydOuLor5EWilUDMG2t74Tbg6iGXZzZ81+q34pc3H7J7iAUA5vGibKrW6gwAuZ8fwsub7Y/eAKxze7V2CxXg6j59vXX25TlSJBc/d28AANxxxx3YunWr+rxz587qvxctWoTFixdj1apVuPXWW/Hiiy/i/vvvx+HDhxESEgIAmD17Nj799FOsW7cOXbt2xZw5c5CamorS0lKrdRE1NzWpN47/67zdP/Rrd1Vh7a4qPHxXNBJ6havLw4ICUPLjWbxbbH9eRPpvemDqsN5O22ZvNDWpN6C5ehjL3gDC37YdxZFT9bjv1ggAV2scExaIlV9X4KOyarvrZJ1tTRjYA8Nu7Yb5nx3EZ+W2dWteZ0uNK/91AZ99W43NB2zn7ABAanw0nk3tx7BDUvCIwOPn52c1qmMhhMDrr7+OZ599FuPGjQMAvPvuu4iKisLatWsxdepUGI1GrFixAu+//z5Gjx4NAFi9ejViYmKwdetWpKSkuLQv5F1mjuyDf+yqcjhv4aOyaocfui1pAMwc9esO2zaZTB3WG4Njw/GQg2vobDlwClscfOi2xDo7Fq0E4tnUfthYXm13n77WOjPskEzcfkgLAI4cOQK9Xo/Y2FhMnDgRP/74IwCgoqICBoMBycnJalutVoukpCQUFl79w1laWoqmpiarNnq9HnFxcWobe8xmM0wmk9WDfE9bh1yuBQ+xtK5/TBiyO6DWrHPrOmqfZp1JNm4PPIMGDcJ7772HzZs3Y/ny5TAYDBgyZAjOnj0Lg8EAAIiKirL6maioKPU1g8GAgIAAhIWFOWxjT25uLhRFUR8xMTEd3DPyFlOTeiP7t33tTq5tr+kjevMQSztMTeqN6SOuv06sc/uwzkS23B54xo4di0ceeQTx8fEYPXo0Nm7cCODqoSsLjcb6o0gIYbOspbbaZGdnw2g0qo+qKvsT9sg3TB3WG4XZI5H+m2u/nsv0Eb3xHykdM0rkC/4jpe91fRizzteGdSay5vbA01JwcDDi4+Nx5MgRdV5Py5GampoaddRHp9OhsbERtbW1DtvYo9VqERoaavUg3xatBGLBuPh2j/Z0ApA9ti8/HK7Df6T0vaZRNdb5+rDORL/wiEnLzZnNZhw8eBD33XcfYmNjodPpkJ+fj7vvvhsA0NjYiB07duCll14CACQkJMDf3x/5+fkYP348AKC6uhr79u3DokWL3NYP8l5Th/VGWn89So/Voq6h0eq1sKAAdA8LxIXGK+gVEcQ5DjfAUZ0tNa76VwM0GuCenmGs8w2wV+fmNa5raERYUADrTNJze+CZO3cuHnzwQfTo0QM1NTV48cUXYTKZMGnSJGg0GsyePRsLFixAnz590KdPHyxYsABBQUFIT08HACiKgilTpmDOnDno2rUrwsPDMXfuXPUQGdH1iFYCkdqff/ydrbU6+/oVlDuSozqzxuRL3B54Tpw4gT/84Q84c+YMunXrhsGDB6O4uBg9e/YEADz11FNoaGjAtGnTUFtbi0GDBmHLli3qNXgA4LXXXoOfnx/Gjx+PhoYGjBo1CqtWreI1eIiIiAgAoBFCOLoEiU8xmUxQFAVGo7HD5vOcbzyPm3JvAgDUZ9cjOCC4Q9brVJfOA/+8us0YXw/4ecE2A9673UREknP2Z2F7P789btIyERERUUdj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUmPgYeIiIikx8BDRERE0mPgISIiIukx8BAREZH0GHiIiIhIegw8REREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9Bh4iIiKSHgMPERERSU+qwPPmm28iNjYWXbp0QUJCAr766it3bxIRERF5AGkCz/r16zF79mw8++yz2Lt3L+677z6MHTsWlZWV7t40IiIicjM/d29AR1m8eDGmTJmCxx57DADw+uuvY/PmzVi6dClyc3PdvHVeRIhf/n3pvPu241p507YSEZHLSRF4GhsbUVpainnz5lktT05ORmFhod2fMZvNMJvN6nOTyeTUbfQaly/88u8Po9y3HURERB1IikNaZ86cweXLlxEVZf0BHRUVBYPBYPdncnNzoSiK+oiJiXHFppKzdbsX6Bzk7q0gIiIPI8UIj4VGo7F6LoSwWWaRnZ2NrKws9bnJZGLoAQBtN2Dcqav/7hwEOKifx/LGbSYiIqeTIvBERESgc+fONqM5NTU1NqM+FlqtFlqt1hWb5100GqBLpLu3goiIqENJcUgrICAACQkJyM/Pt1qen5+PIUOGuGmriIiIyFNIMcIDAFlZWcjIyMCAAQOQmJiIt99+G5WVlXjiiSfcvWlERETkZtIEngkTJuDs2bN44YUXUF1djbi4OGzatAk9e/Z096YRERGRm0kTeABg2rRpmDZtmrs3g4iIiDyMFHN4iIiIiFrDwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUmPgYeIiIikx8BDRERE0mPgISIiIukx8BAREZH0GHiIiIhIegw8REREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikp5bA0+vXr2g0WisHvPmzbNqU1lZiQcffBDBwcGIiIjArFmz0NjYaNWmvLwcSUlJCAwMxC233IIXXngBQghXdoWIiIg8mJ+7N+CFF15AZmam+vymm25S/3358mU88MAD6NatGwoKCnD27FlMmjQJQgi88cYbAACTyYT7778fI0aMQElJCb7//ntMnjwZwcHBmDNnjsv7Q0RERJ7H7YEnJCQEOp3O7mtbtmzBgQMHUFVVBb1eDwB49dVXMXnyZMyfPx+hoaFYs2YNLl68iFWrVkGr1SIuLg7ff/89Fi9ejKysLGg0Gld2h4iIiDyQ2+fwvPTSS+jatSvuuusuzJ8/3+pwVVFREeLi4tSwAwApKSkwm80oLS1V2yQlJUGr1Vq1OXnyJI4dO+bw95rNZphMJqsHERERycmtIzxPPvkk7rnnHoSFhWH37t3Izs5GRUUF/v73vwMADAYDoqKirH4mLCwMAQEBMBgMaptevXpZtbH8jMFgQGxsrN3fnZubi+eff76De0RERESeqMNHeHJycmwmIrd87NmzBwDwl7/8BUlJSbjzzjvx2GOP4a233sKKFStw9uxZdX32DkkJIayWt2xjmbDc2uGs7OxsGI1G9VFVVXVD/SYiIiLP1eEjPDNmzMDEiRNbbdNyRMZi8ODBAIAffvgBXbt2hU6nw65du6za1NbWoqmpSR3F0el06miPRU1NDQDYjA41p9VqrQ6DERERkbw6PPBEREQgIiLiun527969AIDo6GgAQGJiIubPn4/q6mp12ZYtW6DVapGQkKC2eeaZZ9DY2IiAgAC1jV6vdxisiIiIyLe4bdJyUVERXnvtNZSVlaGiogL//Oc/MXXqVKSlpaFHjx4AgOTkZNx+++3IyMjA3r178cUXX2Du3LnIzMxEaGgoACA9PR1arRaTJ0/Gvn37sGHDBixYsIBnaBEREZHKbZOWtVot1q9fj+effx5msxk9e/ZEZmYmnnrqKbVN586dsXHjRkybNg333nsvAgMDkZ6ejldeeUVtoygK8vPzMX36dAwYMABhYWHIyspCVlaWO7pFREREHshtgeeee+5BcXFxm+169OiBzz77rNU28fHx2LlzZ0dtGhEREUnG7dfhISIiInI2Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUmPgYeIiIikx8BDRERE0mPgISIiIukx8BAREZH0GHiIiIhIegw8REREJD0GHiIiIpIeAw8RERFJj4GHiIiIpMfAQ0RERNJj4CEiIiLpMfAQERGR9Bh4iIiISHoMPERERCQ9Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfScGnjmz5+PIUOGICgoCDfffLPdNpWVlXjwwQcRHByMiIgIzJo1C42NjVZtysvLkZSUhMDAQNxyyy144YUXIISwarNjxw4kJCSgS5cu+NWvfoW33nrLWd0iIiIiL+PnzJU3Njbi97//PRITE7FixQqb1y9fvowHHngA3bp1Q0FBAc6ePYtJkyZBCIE33ngDAGAymXD//fdjxIgRKCkpwffff4/JkycjODgYc+bMAQBUVFTgt7/9LTIzM7F69Wp8/fXXmDZtGrp164ZHHnnEmV0kIiIiL+DUwPP8888DAFatWmX39S1btuDAgQOoqqqCXq8HALz66quYPHky5s+fj9DQUKxZswYXL17EqlWroNVqERcXh++//x6LFy9GVlYWNBoN3nrrLfTo0QOvv/46AKBfv37Ys2cPXnnlFbcGHoFfRqHON51323YQERG5i6d8/jk18LSlqKgIcXFxatgBgJSUFJjNZpSWlmLEiBEoKipCUlIStFqtVZvs7GwcO3YMsbGxKCoqQnJystW6U1JSsGLFCjQ1NcHf39/md5vNZpjNZvW5yWTq8P5daLqg/jvqlagOXz8RERG1j1snLRsMBkRFWQeBsLAwBAQEwGAwOGxjed5Wm0uXLuHMmTN2f3dubi4URVEfMTExHdInIiIisnVvzL0I8g9y2++/5hGenJwc9VCVIyUlJRgwYEC71qfRaGyWCSGslrdsY5mwfK1tmsvOzkZWVpb63GQydXjo6RbUDafmngIABPkHQQP720JERCS7IP8gh5/JrnDNgWfGjBmYOHFiq2169erVrnXpdDrs2rXLalltbS2amprUERudTqeO5FjU1NQAQJtt/Pz80LVrV7u/W6vVWh0mcwaNRoPI4Ein/g4iIiJq2zUHnoiICERERHTIL09MTMT8+fNRXV2N6OhoAFcnMmu1WiQkJKhtnnnmGTQ2NiIgIEBto9fr1WCVmJiITz/91GrdW7ZswYABA+zO3yEiIiLf4tQ5PJWVlSgrK0NlZSUuX76MsrIylJWVob6+HgCQnJyM22+/HRkZGdi7dy+++OILzJ07F5mZmQgNDQUApKenQ6vVYvLkydi3bx82bNiABQsWqGdoAcATTzyB48ePIysrCwcPHsQ777yDFStWYO7cuc7sHhEREXkL4USTJk0SAGwe27ZtU9scP35cPPDAAyIwMFCEh4eLGTNmiIsXL1qt57vvvhP33Xef0Gq1QqfTiZycHHHlyhWrNtu3bxd33323CAgIEL169RJLly69pm01Go0CgDAajdfdXyIiInKt9n5+a4RoccliH2UymaAoCoxGozq6RERERJ6tvZ/fvJcWERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUnvmm8eKivLBadNJpObt4SIiIjay/K53daNIxh4fnbu3DkAQExMjJu3hIiIiK7VuXPnoCiKw9d5L62fXblyBSdPnkRISIh6F3ZHTCYTYmJiUFVVJeV9t2TvHyB/H2XvHyB/H9k/7yd7Hz2lf0IInDt3Dnq9Hp06OZ6pwxGen3Xq1Andu3e/pp8JDQ2Vcie2kL1/gPx9lL1/gPx9ZP+8n+x99IT+tTayY8FJy0RERCQ9Bh4iIiKSHgPPddBqtXjuueeg1WrdvSlOIXv/APn7KHv/APn7yP55P9n76G3946RlIiIikh5HeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUmPgednx44dw5QpUxAbG4vAwED07t0bzz33HBobG63aaTQam8dbb71l1aa8vBxJSUkIDAzELbfcghdeeMHmHh87duxAQkICunTpgl/96lc263CnN998E7GxsejSpQsSEhLw1VdfuXuTbOTm5mLgwIEICQlBZGQkHn74YRw+fNiqzeTJk23eq8GDB1u1MZvNmDlzJiIiIhAcHIy0tDScOHHCqk1tbS0yMjKgKAoURUFGRgbq6uqc3UXk5OTYbL9Op1NfF0IgJycHer0egYGBGD58OPbv3+81/evVq5fd/0/Tp08H4H3v386dO/Hggw9Cr9dDo9Hgo48+snrdle9XZWUlHnzwQQQHByMiIgKzZs2y+VvW0X1samrC008/jfj4eAQHB0Ov1+Pf//3fcfLkSat1DB8+3OZ9nThxokf0sa330JX7pDv6Z+//o0ajwcsvv6y28eT3r02ChBBCfP7552Ly5Mli8+bN4ujRo+Ljjz8WkZGRYs6cOVbtAIiVK1eK6upq9XHhwgX1daPRKKKiosTEiRNFeXm5+OCDD0RISIh45ZVX1DY//vijCAoKEk8++aQ4cOCAWL58ufD39xf/93//57L+OrJu3Trh7+8vli9fLg4cOCCefPJJERwcLI4fP+7uTbOSkpIiVq5cKfbt2yfKysrEAw88IHr06CHq6+vVNpMmTRJjxoyxeq/Onj1rtZ4nnnhC3HLLLSI/P1988803YsSIEaJ///7i0qVLapsxY8aIuLg4UVhYKAoLC0VcXJxITU11eh+fe+45cccdd1htf01Njfr6woULRUhIiPjggw9EeXm5mDBhgoiOjhYmk8kr+ldTU2PVt/z8fAFAbNu2TQjhfe/fpk2bxLPPPis++OADAUBs2LDB6nVXvV+XLl0ScXFxYsSIEeKbb74R+fn5Qq/XixkzZji1j3V1dWL06NFi/fr14tChQ6KoqEgMGjRIJCQkWK0jKSlJZGZmWr2vdXV1Vm3c1ce23kNX7ZPu6l/zflVXV4t33nlHaDQacfToUbWNJ79/bWHgacWiRYtEbGys1TJ7O0lzb775plAURVy8eFFdlpubK/R6vbhy5YoQQoinnnpK9O3b1+rnpk6dKgYPHtxxG3+dfvOb34gnnnjCalnfvn3FvHnz3LRF7VNTUyMAiB07dqjLJk2aJB566CGHP1NXVyf8/f3FunXr1GU//fST6NSpk8jLyxNCCHHgwAEBQBQXF6ttioqKBABx6NChju9IM88995zo37+/3deuXLkidDqdWLhwobrs4sWLQlEU8dZbbwkhPL9/LT355JOid+/e6v8Tb37/Wv6dcOX7tWnTJtGpUyfx008/qW3+8Y9/CK1WK4xGo9P6aM/u3bsFAKsvTElJSeLJJ590+DOe0kdHgccV+6S7+tfSQw89JEaOHGm1zFveP3t4SKsVRqMR4eHhNstnzJiBiIgIDBw4EG+99RauXLmivlZUVISkpCSrCzGlpKTg5MmTOHbsmNomOTnZap0pKSnYs2cPmpqanNOZdmhsbERpaanNtiUnJ6OwsNBNW9U+RqMRAGzer+3btyMyMhK33norMjMzUVNTo75WWlqKpqYmq/7q9XrExcWp/S0qKoKiKBg0aJDaZvDgwVAUxSU1OXLkCPR6PWJjYzFx4kT8+OOPAICKigoYDAarbddqtUhKSlK3yxv6Z9HY2IjVq1fjT3/6k9XNe739/bNw5ftVVFSEuLg46PV6tU1KSgrMZjNKS0ud2s+WjEYjNBoNbr75Zqvla9asQUREBO644w7MnTsX586dU1/z9D66Yp/0hPfw1KlT2LhxI6ZMmWLzmre+f7x5qANHjx7FG2+8gVdffdVq+X//939j1KhRCAwMxBdffIE5c+bgzJkz+M///E8AgMFgQK9evax+JioqSn0tNjYWBoNBXda8zaVLl3DmzBlER0c7r2OtOHPmDC5fvmx32wwGg1u2qT2EEMjKysLQoUMRFxenLh87dix+//vfo2fPnqioqMB//dd/YeTIkSgtLYVWq4XBYEBAQADCwsKs1te8vwaDAZGRkTa/MzIy0uk1GTRoEN577z3ceuutOHXqFF588UUMGTIE+/fvV3+3vffq+PHj6rZ7cv+a++ijj1BXV4fJkyery7z9/WvOle+Xvb8vYWFhCAgIcGmfL168iHnz5iE9Pd3qxpKPPvooYmNjodPpsG/fPmRnZ+Pbb79Ffn6+uv2e2kdX7ZOe8B6+++67CAkJwbhx46yWe/P7J33gycnJwfPPP99qm5KSEgwYMEB9fvLkSYwZMwa///3v8dhjj1m1tQQbALjrrrsAAC+88ILV8ubfUAGoE5abL29PG3ext22esF2OzJgxA9999x0KCgqslk+YMEH9d1xcHAYMGICePXti48aNNv+Jm2vZX3t9d0VNxo4dq/47Pj4eiYmJ6N27N9599111ouT1vFee0r/mVqxYgbFjx1p94/P2988eV71f7u5zU1MTJk6ciCtXruDNN9+0ei0zM1P9d1xcHPr06YMBAwbgm2++wT333APAc/voyn3S3e/hO++8g0cffRRdunSxWu7N75/0h7RmzJiBgwcPtvpoPipw8uRJjBgxAomJiXj77bfbXP/gwYNhMplw6tQpAIBOp7NJqJYhT0uiddTGz88PXbt2vaH+3oiIiAh07tzZ7ra1TOOeYubMmfjkk0+wbds2dO/evdW20dHR6NmzJ44cOQLg6vvQ2NiI2tpaq3bN+6vT6dT3trnTp0+7vCbBwcGIj4/HkSNH1LO1WnuvvKV/x48fx9atW22+XLTkze+fK98ve39famtr0dTU5JI+NzU1Yfz48aioqEB+fr7V6I4999xzD/z9/a3eV0/vo4Wz9kl39++rr77C4cOH2/w/CXjZ++e02UFe6MSJE6JPnz5i4sSJVjPqW/PGG2+ILl26qJOU33zzTXHzzTcLs9mstlm4cKHNpOV+/fpZreeJJ57wmEnLf/7zn62W9evXz+MmLV+5ckVMnz5d6PV68f3337frZ86cOSO0Wq149913hRC/TDBcv3692ubkyZN2Jxju2rVLbVNcXOyWSb0XL14Ut9xyi3j++efVSbAvvfSS+rrZbLY7CdbT+/fcc88JnU4nmpqaWm3nTe8fHExadsX7ZZkQevLkSbXNunXrXDJpubGxUTz88MPijjvusDqjsDXl5eVWJxx4Sh/t9a8lZ+2T7u7fpEmTbM6uc8RT3z97GHh+9tNPP4lf//rXYuTIkeLEiRNWp9xZfPLJJ+Ltt98W5eXl4ocffhDLly8XoaGhYtasWWqburo6ERUVJf7whz+I8vJy8eGHH4rQ0FC7p6X/5S9/EQcOHBArVqzwuNPSV6xYIQ4cOCBmz54tgoODxbFjx9y9aVb+/Oc/C0VRxPbt2+1eIuDcuXNizpw5orCwUFRUVIht27aJxMREccstt9icBty9e3exdetW8c0334iRI0faPYX0zjvvFEVFRaKoqEjEx8e75LTtOXPmiO3bt4sff/xRFBcXi9TUVBESEqK+FwsXLhSKoogPP/xQlJeXiz/84Q92T3P21P4JIcTly5dFjx49xNNPP2213Bvfv3Pnzom9e/eKvXv3CgBi8eLFYu/eveoZSq56vyyn/I4aNUp88803YuvWraJ79+4dcspva31samoSaWlponv37qKsrMzq/6XlC+APP/wgnn/+eVFSUiIqKirExo0bRd++fcXdd9/tEX1srX+u3Cfd0T8Lo9EogoKCxNKlS21+3tPfv7Yw8Pxs5cqVAoDdh8Xnn38u7rrrLnHTTTeJoKAgERcXJ15//XWbb6bfffeduO+++4RWqxU6nU7k5OSoozsW27dvF3fffbcICAgQvXr1srtzucvf/vY30bNnTxEQECDuueceq1O9PYWj92rlypVCCCEuXLggkpOTRbdu3YS/v7/o0aOHmDRpkqisrLRaT0NDg5gxY4YIDw8XgYGBIjU11abN2bNnxaOPPipCQkJESEiIePTRR0Vtba3T+2i5Tou/v7/Q6/Vi3LhxYv/+/errV65cUUdHtFqtGDZsmCgvL/ea/gkhxObNmwUAcfjwYavl3vj+bdu2ze4+OWnSJCGEa9+v48ePiwceeEAEBgaK8PBwMWPGDKtLZTijjxUVFQ7/X1qurVRZWSmGDRsmwsPDRUBAgOjdu7eYNWuWzbVs3NXH1vrn6n3S1f2zWLZsmQgMDLS5to4Qnv/+tUUjRItLABMRERFJRvpJy0REREQMPERERCQ9Bh4iIiKSHgMPERERSY+Bh4iIiKTHwENERETSY+AhIiIi6THwEBERkfQYeIiIiEh6DDxEREQkPQYeIiIikh4DDxEREUnv/wOPrwuB2YfcBgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[1000])\n",
    "\n",
    "output_file = pset.ParticleFile(\n",
    "    name=\"NestedFieldParticle.zarr\", outputdt=50, chunks=(1, 100)\n",
    ")\n",
    "\n",
    "pset.execute(AdvectionRK4, runtime=14000, dt=10, output_file=output_file)\n",
    "\n",
    "ds = xr.open_zarr(\"NestedFieldParticle.zarr\")\n",
    "plt.plot(ds.lon.T, ds.lat.T, \".-\")\n",
    "plt.plot([0, 2e3, 2e3, 0, 0], [0, 0, 2e3, 2e3, 0], c=\"orange\")\n",
    "plt.plot([-2e3, 18e3, 18e3, -2e3, -2e3], [-1e3, -1e3, 3e3, 3e3, -1e3], c=\"green\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we observe, there is a change of dynamic at lon=2000, which corresponds to the change of grid.\n",
    "\n",
    "The analytical solution to the problem:\n",
    "\n",
    "\\begin{align}\n",
    "dx/dt &= 1;\\\\\n",
    "dy/dt &= \\cos(x \\pi/400);\\\\\n",
    "\\text{with } x(0) &= 0, y(0) = 1000\n",
    "\\end{align}\n",
    "\n",
    "is\n",
    "\n",
    "\\begin{align}\n",
    "x(t) &= t;\\\\\n",
    "y(t) &= 1000 + 400/\\pi \\sin(t \\pi / 400)\n",
    "\\end{align}\n",
    "\n",
    "which is captured by the High Resolution field (orange area) but not the Low Resolution one (green area).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Keep track of the field interpolated\n",
    "\n",
    "For different reasons, you may want to keep track of the field you have interpolated. You can do that easily by creating another field that share the grid with original fields.\n",
    "Watch out that this operation has a cost of a full interpolation operation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need to redefine fieldset\n",
    "fieldset = FieldSet(U, V)\n",
    "\n",
    "ones_array1 = np.ones((U1.grid.ydim, U1.grid.xdim), dtype=np.float32)\n",
    "F1 = Field(\"F1\", ones_array1, grid=U1.grid)\n",
    "\n",
    "ones_array2 = np.ones((U2.grid.ydim, U2.grid.xdim), dtype=np.float32)\n",
    "F2 = Field(\"F2\", 2 * ones_array2, grid=U2.grid)\n",
    "\n",
    "F = NestedField(\"F\", [F1, F2])\n",
    "fieldset.add_field(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1.0/1.0 [00:00<00:00, 528.52it/s]\n",
      "Particle (1000, 500) interpolates Field #1\n",
      "100%|██████████| 1.0/1.0 [00:00<00:00, 1368.01it/s]\n",
      "Particle (1000, 500) interpolates Field #1\n"
     ]
    }
   ],
   "source": [
    "def SampleNestedFieldIndex(particle, fieldset, time):\n",
    "    particle.f = fieldset.F[time, particle.depth, particle.lat, particle.lon]\n",
    "\n",
    "\n",
    "SampleParticle = JITParticle.add_variable(\"f\", dtype=np.int32)\n",
    "\n",
    "pset = ParticleSet(fieldset, pclass=SampleParticle, lon=[1000], lat=[500])\n",
    "pset.execute(SampleNestedFieldIndex, runtime=1)\n",
    "print(\n",
    "    f\"Particle ({pset[0].lon:g}, {pset[0].lat:g}) \"\n",
    "    f\"interpolates Field #{int(pset[0].f)}\"\n",
    ")\n",
    "\n",
    "pset[0].lon = 10000\n",
    "pset.execute(SampleNestedFieldIndex, runtime=1)\n",
    "print(\n",
    "    f\"Particle ({pset[0].lon:g}, {pset[0].lat:g}) \"\n",
    "    f\"interpolates Field #{int(pset[0].f)}\"\n",
    ")"
   ]
  }
 ],
 "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}