{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## MITgcm Example \n", "\n", "xgcm is developed in close coordination with the [xmitgcm](http://xmitgcm.readthedocs.io/) package.\n", "The metadata in datasets constructed by xmitgcm should always be compatible with xgcm's expectations.\n", "xmitgcm is necessary for reading MITgcm's binary MDS file format.\n", "However, for this example, the MDS files have already been converted and saved as netCDF.\n", "\n", "Below are some example of how to make calculations on mitgcm-style datasets using xgcm.\n", "\n", "First we import xarray and xgcm:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#! pip install git+https://github.com/xgcm/xgcm.git" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import numpy as np\n", "import xgcm\n", "from matplotlib import pyplot as plt\n", "\n", "%matplotlib inline\n", "plt.rcParams['figure.figsize'] = (10,6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we open the example dataset, which is stored with the xgcm github repository in the datasets folder." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "#!wget http://www.ldeo.columbia.edu/~rpa/mitgcm_example_dataset.nc" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['__class__',\n", " '__delattr__',\n", " '__dict__',\n", " '__dir__',\n", " '__doc__',\n", " '__eq__',\n", " '__format__',\n", " '__ge__',\n", " '__getattribute__',\n", " '__gt__',\n", " '__hash__',\n", " '__init__',\n", " '__init_subclass__',\n", " '__le__',\n", " '__lt__',\n", " '__module__',\n", " '__ne__',\n", " '__new__',\n", " '__reduce__',\n", " '__reduce_ex__',\n", " '__repr__',\n", " '__setattr__',\n", " '__sizeof__',\n", " '__str__',\n", " '__subclasshook__',\n", " '__weakref__',\n", " '_apply_vector_function',\n", " '_assign_face_connections',\n", " 'cumsum',\n", " 'diff',\n", " 'diff_2d_vector',\n", " 'interp',\n", " 'interp_2d_vector']" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dir(xgcm.grid.Grid)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "\n", "Dimensions: (XC: 90, XG: 90, YC: 40, YG: 40, Z: 15, Zl: 15, Zp1: 16, Zu: 15, layer_1RHO_bounds: 31, layer_1RHO_center: 30, layer_1RHO_interface: 29, time: 1)\n", "Coordinates:\n", " iter (time) int64 ...\n", " * time (time) timedelta64[ns] 11:00:00\n", " * XC (XC) float32 2.0 6.0 10.0 14.0 18.0 22.0 26.0 30.0 ...\n", " * YC (YC) float32 -78.0 -74.0 -70.0 -66.0 -62.0 -58.0 ...\n", " * XG (XG) float32 0.0 4.0 8.0 12.0 16.0 20.0 24.0 28.0 ...\n", " * YG (YG) float32 -80.0 -76.0 -72.0 -68.0 -64.0 -60.0 ...\n", " * Z (Z) float32 -25.0 -85.0 -170.0 -290.0 -455.0 ...\n", " * Zp1 (Zp1) float32 0.0 -50.0 -120.0 -220.0 -360.0 ...\n", " * Zu (Zu) float32 -50.0 -120.0 -220.0 -360.0 -550.0 ...\n", " * Zl (Zl) float32 0.0 -50.0 -120.0 -220.0 -360.0 -550.0 ...\n", " rA (YC, XC) float32 ...\n", " dxG (YG, XC) float32 ...\n", " dyG (YC, XG) float32 ...\n", " Depth (YC, XC) float32 ...\n", " rAz (YG, XG) float32 ...\n", " dxC (YC, XG) float32 ...\n", " dyC (YG, XC) float32 ...\n", " rAw (YC, XG) float32 ...\n", " rAs (YG, XC) float32 ...\n", " drC (Zp1) float32 ...\n", " drF (Z) float32 ...\n", " PHrefC (Z) float32 ...\n", " PHrefF (Zp1) float32 ...\n", " hFacC (Z, YC, XC) float32 ...\n", " hFacW (Z, YC, XG) float32 ...\n", " hFacS (Z, YG, XC) float32 ...\n", " * layer_1RHO_bounds (layer_1RHO_bounds) float32 19.9499 20.4499 ...\n", " * layer_1RHO_center (layer_1RHO_center) float32 20.1999 20.6922 ...\n", " * layer_1RHO_interface (layer_1RHO_interface) float32 20.4499 20.934498 ...\n", "Data variables:\n", " Ttave (time, Z, YC, XC) float32 ...\n", " UStave (time, Z, YC, XG) float32 ...\n", " wVeltave (time, Zl, YC, XC) float32 ...\n", " VTtave (time, Z, YG, XC) float32 ...\n", " VVtave (time, Z, YG, XC) float32 ...\n", " V (time, Z, YG, XC) float32 ...\n", " VStave (time, Z, YG, XC) float32 ...\n", " PHL (time, YC, XC) float32 ...\n", " TFLUX (time, YC, XC) float32 ...\n", " surForcT (time, YC, XC) float32 ...\n", " SFLUX (time, YC, XC) float32 ...\n", " surForcS (time, YC, XC) float32 ...\n", " vVeltave (time, Z, YG, XC) float32 ...\n", " GM_Kwz-T (time, Zl, YC, XC) float32 ...\n", " W (time, Zl, YC, XC) float32 ...\n", " uVeltave (time, Z, YC, XG) float32 ...\n", " TOTSTEND (time, Z, YC, XC) float32 ...\n", " ADVr_SLT (time, Zl, YC, XC) float32 ...\n", " ADVx_SLT (time, Z, YC, XG) float32 ...\n", " ADVy_SLT (time, Z, YG, XC) float32 ...\n", " DFrE_SLT (time, Zl, YC, XC) float32 ...\n", " DFxE_SLT (time, Z, YC, XG) float32 ...\n", " DFyE_SLT (time, Z, YG, XC) float32 ...\n", " DFrI_SLT (time, Zl, YC, XC) float32 ...\n", " USLTMASS (time, Z, YC, XG) float32 ...\n", " VSLTMASS (time, Z, YG, XC) float32 ...\n", " WSLTMASS (time, Zl, YC, XC) float32 ...\n", " tFluxtave (time, YC, XC) float32 ...\n", " LaTs1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LaTh1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LaTz1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LTha1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LTza1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LaSs1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LaSh1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LaSz1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LSha1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LSza1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LTto1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " LSto1RHO (time, layer_1RHO_interface, YC, XC) float32 ...\n", " WTtave (time, Zl, YC, XC) float32 ...\n", " LaUH1RHO (time, layer_1RHO_center, YC, XG) float32 ...\n", " LaHw1RHO (time, layer_1RHO_center, YC, XG) float32 ...\n", " LaPw1RHO (time, layer_1RHO_center, YC, XG) float32 ...\n", " LaVH1RHO (time, layer_1RHO_center, YG, XC) float32 ...\n", " LaHs1RHO (time, layer_1RHO_center, YG, XC) float32 ...\n", " LaPs1RHO (time, layer_1RHO_center, YG, XC) float32 ...\n", " U (time, Z, YC, XG) float32 ...\n", " Convtave (time, Zl, YC, XC) float32 ...\n", " sFluxtave (time, YC, XC) float32 ...\n", " PHLtave (time, YC, XC) float32 ...\n", " GM_Kwy-T (time, Zl, YC, XC) float32 ...\n", " uFluxtave (time, YC, XG) float32 ...\n", " GM_Kwx-T (time, Zl, YC, XC) float32 ...\n", " WStave (time, Zl, YC, XC) float32 ...\n", " S (time, Z, YC, XC) float32 ...\n", " UUtave (time, Z, YC, XG) float32 ...\n", " TTtave (time, Z, YC, XC) float32 ...\n", " UTtave (time, Z, YC, XG) float32 ...\n", " Eta2tave (time, YC, XC) float32 ...\n", " Stave (time, Z, YC, XC) float32 ...\n", " PhHytave (time, Z, YC, XC) float32 ...\n", " ETAtave (time, YC, XC) float32 ...\n", " UVEL (time, Z, YC, XG) float32 ...\n", " VVEL (time, Z, YG, XC) float32 ...\n", " WVEL (time, Zl, YC, XC) float32 ...\n", " THETA (time, Z, YC, XC) float32 ...\n", " SALT (time, Z, YC, XC) float32 ...\n", " PHL2tave (time, YC, XC) float32 ...\n", " T (time, Z, YC, XC) float32 ...\n", " Eta (time, YC, XC) float32 ...\n", " UVtave (time, Z, YG, XG) float32 ...\n", " Tdiftave (time, Zl, YC, XC) float32 ...\n", " vFluxtave (time, YG, XC) float32 ...\n", " PH (time, Z, YC, XC) float32 ...\n", " TOTTTEND (time, Z, YC, XC) float32 ...\n", " ADVr_TH (time, Zl, YC, XC) float32 ...\n", " ADVx_TH (time, Z, YC, XG) float32 ...\n", " ADVy_TH (time, Z, YG, XC) float32 ...\n", " DFrE_TH (time, Zl, YC, XC) float32 ...\n", " DFxE_TH (time, Z, YC, XG) float32 ...\n", " DFyE_TH (time, Z, YG, XC) float32 ...\n", " DFrI_TH (time, Zl, YC, XC) float32 ...\n", " UTHMASS (time, Z, YC, XG) float32 ...\n", " VTHMASS (time, Z, YG, XC) float32 ...\n", " WTHMASS (time, Zl, YC, XC) float32 ...\n", "Attributes:\n", " Conventions: CF-1.6\n", " title: netCDF wrapper of MITgcm MDS binary data\n", " source: MITgcm\n", " history: Created by calling `open_mdsdataset(llc_method='smallchunks..." ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds = xr.open_dataset('mitgcm_example_dataset.nc')\n", "ds" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "array([[[ 0. , 0. , ..., 0. , 0. ],\n", " [ 0. , 0. , ..., 0. , 0. ],\n", " ...,\n", " [-0.702411, -0.694586, ..., -0.694789, -0.709546],\n", " [-0.656158, -0.659816, ..., -0.623395, -0.649 ]]], dtype=float32)\n", "Coordinates:\n", " iter (time) int64 ...\n", " * time (time) timedelta64[ns] 11:00:00\n", " * XC (XC) float32 2.0 6.0 10.0 14.0 18.0 22.0 26.0 30.0 34.0 38.0 ...\n", " * YC (YC) float32 -78.0 -74.0 -70.0 -66.0 -62.0 -58.0 -54.0 -50.0 ...\n", " rA (YC, XC) float32 ...\n", " Depth (YC, XC) float32 ...\n", "Attributes:\n", " standard_name: ETAN\n", " long_name: Surface Height Anomaly\n", " units: m" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds.Eta" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXucH3V579/Pbi5LkiVhExKSQFwQPBQQA0b01GpRUOOlpR61B2/FFkutl9Z6rELtUatioRepx6qYKor1gohaaOsNUbS2VQgakEiUgGtIWBOTSNgEQkj2OX/MzGZ2du733/ye9+v1e+38Zr4z8/xmZ76f7/N8n+93RFUxDMMwjCIMNG2AYRiG0fuYmBiGYRiFMTExDMMwCmNiYhiGYRTGxMQwDMMojImJYRiGURgTkxhEZKOInN20HV1FRF4uIl9v2o6qEJF3isinmrbDMOrAxCQGVT1VVW+GdlUMIvJGEblXRB4UkftF5AoRmeXb/usicouITIjIHSLyG4H9jxaRz4jIAyLyKxH5tG/bXBG5yj32L0TkTYF9zxGRTSLykIh8S0Qek9LmURFRv52q+mlVfXb+K1E9IrJaRG5zf+9tIrK6aZuiEJE5InKdiIy51/rswPZnuP+zPSIyluJ4o275h9z/+bmB7X/m3iN73Htmbh3HMtqJiUlN+CvREvhX4ExVPRI4DXgC8CfueUaAG4C/BRYBfwP8q4gc5dv/i8AvgMcAS4G/8217J3CSu+0ZwFtEZK177CXuvv8XGAHWA58r8Xe1ChGZA1wPfAo4CrgauN5d31a+C7wC5/8bZB9wFfDnKY/1WeCHwGLgbcB1InI0gIg8B7gYOAcYBU4A/qqmYxltRFXtE/EBxoBzgbXAAeBRYC9wu7t9IfAxYBzYBrwHGHS3vQr4T+AKYDfwnopsXAx8A/iQ+/0FwMZAmZ8CF7rLz3Z/12DE8bYBz/Z9fzdwjbt8EfBfvm3zgYeBk1PYuQVQ9/rtBf6ne42+6yujwGuBu4EJ99yPBf4beBC4FpjjK/8CYAPwAPBfwOklX9tnu9dDAr9jbUT544Fvu7bfCPwj8Cnf9s/jVPJ7gO8Ap7rrnwRsB2b5yr4I2OAun4Uj3A+65d6XwvatwNkR284FxhL2fxzwCDDsW/cfwGvc5c8A7/VtOwf4RdXHsk97P+aZpEBVvwq8F/icqi5Q1Se4m64GDgInAmfgVD6v9u36ZOBenNb/pcHjisjL3FBT1GdVlE3uvg8CO3E8k494m9zPtOI4HgzAU4CfAFeLyC4RuVVEftM95lHACuB23763A6e6y6f6t6nqPuAe3/Y4nu7+XeRew/+OKLcWeKJr51uAdcDLgePc3/BS19YzcVrZf4QjqB8BbogKj7jhvqjr/KEIW04F7lC3hnO5I+b3fga4DViCI4QXBLZ/BcfrWwr8APg0gKreCuwCnuUr+wrgn93l9wPvV8cTfSyOqJaOiPybiFzsfj0VuFdVJ3xFIu8Fd3mZiCwu+1hGb1Bm6KWvEJFlwHNxKseHgX0icgVO692r2O9X1Q+4yweDx1DVz+BUQJnx9hWRk4Dfw2mxgtNCXyEiLwWuA16GUwHNc7cfy2HR+32cFvD1InIicIRbZo/vVHuAYXd5AfDLgCn+7WVwuao+CGwUkTuBr6vqvQAi8hUc0b4a+EPgI6r6fXe/q0XkL3BE6NvBg6rq6TlsWcD0awERv9cV/icB56rqI8B3RORfAzZc5Sv/TuBXIrJQVfe4v+kVwFfcUOVzcLw0cDziE0VkiaruBL6X47ckoqov8H2N+u0rI7Z7y8PArjKPldZ+o1nMM8nPY4DZwLjXwsURkaW+MvdVbYSq3g1sBD7kft8FnAe8CUdg1uKEwba6uzyME+L4mKo+qqrXuHY+FSf8BHCk7xRH4oRtcLf7twW3l8F23/LDId8XuMuPAf6P38PA8V5WlGhLlt+7AviV6615/NxbEJFBEblMRO5xPcoxd9MS9++ngN8SkQXA7wL/oarj7rYLcUJFm1xP0l9RV0XSbw9u95bDrk2ZxzJaiolJeoLTK9+HEwdeoqqL3M+RqnpqzD7TECc1dm/MJzLMFWAWjvfhnFT126r6JFUdAV4J/A/gFnfzHVF2qeqvcPp/nuBb/QQcscL9O7VNROa7591IMmVPT30fcKnv2i9S1Xmq+tmwwuKkeUdd5ysjzrEROF1E/GHD0wn/vePAUe418fD//16GI/Ln4vS1jXqmAajqNpy+oRfi/M+8EBeqereqvhSnoXI5Tue1/zxVsBE4QUT8XljkveAub3cbM1Uey2gpJibp2Q6MisgAgNtq/Drw9yJypIgMiMhjvf6HNKiTGrsg5rMlbD8RebWILHWXTwEuAW7ybT9DRGaLyJE4mVpbVfVr7uYv4VR6F7it5RfjhBv+093+SeAvReQoETkZJ5z0Cd++p4nIi0RkCHg7Tp/CJve87xSRmyN+7i+BSZxMnTL4J+A1IvJkcZgvIs8PVFhTqJPmHXWdXxNxjpuBQ8CfiJMy/Xp3/TdDjv9znE7yvxInRfc3gN/yFRnGaXzswgk5vjfkfJ/E6Sd6PM61BkBEXiEiR6vqJE6yAa5dM3DtHHK/zhGRIU8M3Xt0CMejFndbaGaaqv4UJ7nhHW65F+II6Rd8tl4oIqe4fW1/yeH7pLJjGS2m6QyANn9ws7nc5cU4aZe/An7grlsIfBgnhLQHJ/XxfHfbq/BlKpVs18dxxG2fa+PfAkO+7Z917dmDk7q7NLD/04Af4YQX1gNP822bi9Ox7WUOvSmw77nAJpyQ083AqG/bx3C8hSi734UjKg/g9G1Mu0Y43suJvu/fBV7l+/4e4KO+72uBW93jjeNkSw1HnT/ntT4Dp1P9YZxO8zNiyp6Ak6W0l0A2F0547nqc0M3Pcfq5gr93nnvdrw4c91PADve4G4HfSbhnNfAZdbedHbLtZt++XwH+wvd91P0fP4yTtHFu4FxeKPVB956cW8Wx7NMbH3H/kYZRGBHZAJyjFp7IjYjcA/yRqn6jaVsMIwuWzWWUhqq2dnR4LyAiL8LxFmaE0Qyj7ZiYGEYLcPuaTgFeqU7fiGH0FBbmMgzDMApj2VyGYRhGYToR5pq/cESPOmZlckGXObMGmD97kEGflO4/OMmBQ+3y0uYMCrMHBqbZaRiGw6FJeHQy/Lk9cDA6Urjtp3fuVNWji5z7ODlC95MuGrmTA19T1bVFztcLdEJMjjpmJa+/8oupy4+OzGPVwiFOXuzMHrJ7/yF27DvAlj37qzIRgFULhyK3hZ171cIhls6fw8jQYJVmGUZPsnv/Idbf/2DotrHdD0Xud8kzH/fzyI0p2c8kL2J5qrIf4edLkkv1Pp0Qk6x4N9rS+c54raaFxNtetQ2GkZfd+0PHSE6j7kbPjn0HQtfHCYlRHX0pJh479h2YEpQkkm7Q0ZF5sduTzhP1YBhG03ieexIjQ0cklimLNOJm1EvfiolfHMI8gqytm7HdD0UKiueVxLfc5kTaYhhNsmPfAb7zs92J5ZbOP6YW78QTkjKeW6M8+lZMPII3ZJGbMUpQ0vR7ONvnmJgYrWL3/kNs2bOfTePJE/iuH5nHs084KrFcGYR5SiYkzdLXYlLHzZfUV5JU3jrfjSbZse8AY7sfYmx7spiMLR9mky+xpQqiQm4mJM1jSacVkjUba2RocFrfStr+HMNoA5vGJ9iyZ38t/Rl+D96EpB30lWeS1EkOxW7MsOPn8Sw8ETGvxOglxrZPMLZ8mFULhyrpjA9L4TchaQ99ISZpRCRYNkv2VlhfieeVZMUExOhlNo1PMDoyr5LxUX4hMRFpH50VkywCkrS/d+NGHTNqvQmD0W9U5Z34Q2cmJO2kk2JSVEiKHi+vV2IYbWPp/DmMjsxjdNn0F1jGdciX7Z34w1ttEZK5A8Jjj0j5jO+r1pa20EkxCZI1o6ooNgWK0RVGhgadqYeWz3wbcpSgeN7JmhVHlmZHVMp82rC0UT2dE5OwvgvDMKLxh5DCGkGed+Jn0/gEo8uGU6UMZ7UhSNzoe+/5NjFpnk6JSZSQ1B1yMq/E6BWC4zbC+jk878TPycuHSxWUOMGI80q8Z3t0ZJ4JSsN0Skz8+IXEKnfDCCeYahv1vAQbZKMj86ZGxRcVFL9XkmUGCL/ArVo4ZGLSMJ0Rk7I73Q2j63iVuL8SjsrCCgqM148SNc1KlucxbuxIXAbldOGbY95Jw3RGTPyYV2IYyXiVuF8QysrCGtv9UKoOeG/uL78IhAlUlxuLInIV8AJgh6qeFrL9PODdwCRwEHijqn7X3XYI+JFbdIuq/nY9Vs+kc9OpNNVPYhi9RFglDo4IJE03n2byx03jE5EvrvITPFfwmN73sd0PTdm6aXxihp3B44yOzOslAfoEEPcmxpuAJ6jqauAPgI/6tj2sqqvdT2NCAh3xTObMCtdE80oM4zD+vglvAsewytvfsR2Gt69HWH+Jlx4cN/GjX9DSzEocZMue/VN2tmkMSlZU9TsiMhqzfa/v63ygXe8Xd+mEmHhYeMswwglmbcV5FmO7H0p8xXSayj9u8KJnTxoB2DQ+MTXOJXjeJC+qJf0oS0Rkve/7OlVdl+UAIvJC4K+BpcDzfZuG3GMfBC5T1X8pbG1OOiUmhmGEk2WCRE8EokjySvzb4qZW8ezJ45V4djQ1zmRoQDhxwex0hfexU1XXFDmfqn4J+JKIPB2n/+Rcd9MqVb1fRE4AvikiP1LVe4qcKy+NiomILMKJ/52G47r9AfAT4HPAKDAG/K6q/irpWOaVGFWyYDJ7ZbV3oB0x+7Sd3P5R7klik5Yw7yTOK/HEKTh9S9Q5w1KJg0LYEu+kFNyQ2GNFZImq7lTV+93194rIzcAZQCNi0nQH/PuBr6rqycATgLuAi4GbVPUknI6nixu0zzByk0eAqsALBW0an5j6BAmbLiUKf9lgpR9kbHt4Z3mYV5JnrIq/Yz4Or0O+hzrlpxCRE0VE3OUzcd7xvUtEjhKRue76JcBTgR83ZWdjnomIHAk8HXgVgKoeAA64aXBnu8WuBm4G3hp3rDmDUpWZhtEaUagKTxyyVLR5vRNI7ixPEqg4elQsPotT5y0Rka3AO4DZAKp6JfAi4PdE5FHgYeB/q6qKyK8BHxGRSRzH4DJV7T8xAU4Afgl8XESeANwG/CmwTFXHAVR1XESWhu0sIhcBFwEsW3FsPRYbfUdRIVkw+VBrwl1BpnkYI/MyzWMXN2AxjLD+jaJeSVdQ1ZcmbL8cuDxk/X8Bj6/Krqw0KSazgDOBN6jq90Xk/WQIabnZEOsATj59tYKNLTGqY2BiO5PDy5o2ozSCYa2sr03IMp3KVD+Iz2uoQkiivBKb7LUemhSTrcBWVf2++/06HDHZLiLLXa9kObCjMQuNvibolaQVlIGJ7QBTZZv2TpbOn8OWPfsj+0XyjHpPmk4ljDQeSZEQVxgmJPXRWAe8qv4CuE9E/oe76hyczqMbgAvcdRcA1ycda/bAgHklLWf3/kOpPk2zYPKhqY+HJw7B5TCiyra138ULb2XNgAxOS59GBKI6/7McA8KTBcK8kjAhsXqiOpoeZ/IG4NMiMge4F/h9HIG7VkQuBLYAL8lyQEsLbicjQ4OtEIs4oir8yeFlM7yNKLKUrYuwKeT95Klg/S/NShPuGts+keltjUHiss2ShMQEpB4aFRNV3QCEDeY5J+uxTESMqkkrDn5B8dNkuCuuQs377Pi9k6zT0dfV4W5CUh9NeyalMNj0aBmj50kKQ2X1MtrilXg4glFuxRr0eLJ2yKclr1fiCUkVDc1ZcwZZ/JiF6QrHR0Y7QyfExOgNqg51BY/dRm+1Se+kiuvhVdjBcSNhgpLHG4lLGgijDiExwjExMWolTlDyPPhx4pT0bnOPKjvH/eEuf3aXR1vHoKTF83iC4S443KGeRkTCOt+zCkm4bUZdmJgYrWH3/kOZKoA0Xk7bK5Sm04bLwAt3JU1rkiXtt6iQQPb7ySiG9TYYtVPGA16mkNRRmbetD6VsvM740LTdZcO1CUlw4se2ZxB2CRMToxGqbDGODA22pkWaRkTaOgYlC0npx0XJOudW0ntOjPIxMTEaI2+lH9fnkldEsngnWT2ZMEEZmNg+9ekKcd5JGk5ePpx6QGIcnnfiCYp5J/VgYmI0ThleRJWeyN6BedM+RYgSkH73TsroIwnDBKU+rAPeaAVp04brSv8tux+lSx5IGaTxXooIyZY9+6eEbce+AzZ4sQZMTIzWkEYYqhKPoGfgffeLShe8hyrx3uiYRJKQ9MI7SQbnDnLUCYvSFb6lWlvagomJ0ffEiUQZAtIvXon/VbxZZxSG9CISFkpLI2JGtZiYGEZF9IuIQDGvJElE0vTD2FTzzWNiYvQ1VYWu+klIYLpX4hH0ToJvdkzCBKK3sGwuwyiZfhOS3fsP8Z2f7T480ePIvCmx8AQkbbrwqoVDUx+jtzDPxOgr6pqHq59IGiAYFJK6X69rmVz1YGJidJa6sq/6VUQgfV+JR9o3IpZNW2ZE6DImJkajhFX4Wcd4NJWy268iEhzr48/g8n9Pog4R6QWvRETWAu8HBoGPquplge2PAa4CjgZ2A69Q1a3utguAv3SLvkdVr67N8AAmJkbraPN4jn4VEI/d+w/lnvfK75XU2SfSZq9ERAaBDwLPArYCt4rIDar6Y1+xvwM+qapXi8gzgb8GXikiI8A7cN5Wq8Bt7r6/qvdXOJiYGI3RZtHw6HfxCLJj34FcYzrSjiFJ8iQ6OIHjWcBmVb0XQESuAc4D/GJyCvBn7vK3gH9xl58D3Kiqu919bwTWAp+twe4ZNC4mrjKvB7ap6gtE5HjgGmAE+AHwSlXt3B3U77RNSEw0kvH6R7ww1ujIvFQhraCQBL2SLKGopfPnpBaUpfPnVOaVDM6dzcLR1K8VWCIi633f16nqOnd5JXCfb9tW4MmB/W8HXoQTCnshMCwiiyP2XZnWqLJpXEyAPwXuAo50v18OXKGq14jIlcCFwIebMs4oH7+QJFXiRd8DYiJRHjv2HZh6CZYnEFmmPikiIsH9esxD2amqayK2Scg6DXx/M/CPIvIq4DvANuBgyn1ro1ExEZFjgecDlwJvEhEBngm8zC1yNfBOTEw6gyckaSv5nhKDXduity1urMFYCv5O96xzZ5UlIsFj9JigRLEVOM73/Vjgfn8BVb0f+F8AIrIAeJGq7hGRrcDZgX1vrtLYOJr2TP4BeAvgJaIvBh5Q1YPu90i3TUQuAi4COPa448KKGC2lpwQiDXEi4i/TY4KSNItz1k70sjOrOiIotwInueH9bcD5HG5MAyAiS4DdqjoJXIKT2QXwNeC9InKU+/3Z7vZGaExMROQFwA5VvU1EzvZWhxQNddvcmOM6gNVnntmYa2ekZ8HkQ9OFJE0lDO2rhNPa3cN4QhJWWZchIlF9GVnfO9LrgqKqB0Xk9TjCMAhcpaobReRdwHpVvQHH+/hrEVGcMNfr3H13i8i7cQQJ4F1eZ3wTNOmZPBX4bRF5HjCE02fyD8AiEZnleiczXD6jN5nW4Z61Ms5avmzxKUs8esw7Saqks3oaWV8x0C8vtFLVLwNfDqx7u2/5OuC6iH2v4rCn0iiNiYmqXoLrkrmeyZtV9eUi8nngxTgZXRcA1zdlo5GduCytuPDWwfGxqeVZy0eLGdFmzyFGUBZMPlT6S7ny4I0l8acA+72RKkTE6H2a7jMJ463ANSLyHuCHwMcatsdISap0X19F7xcQP6UKi1EYT0iCImIiYfhphZio6s24WQju4J2zmrTHKIcoTyRKRLKU7WmRCXgnAxPbp1Kg2+KdeASFxATEiKIVYmL0PqnGjuzaFikOj2zbAsDclatSnS9OkHpaaFpIlJAEPdE2iaBRPyYmRmHSCokfTzyC+NenFZYgYULTSwLTFu9k1cKhGaPIy3zFcRt+Y14G585m4Ym9k0xRByYmRi7CKo6ktF+vko8SkiDBcnnFxX9uP20SGH+oq62UPQVOkmiODA0mZnT1clpw1zAxMTJTh5CEkTUUlkSjHf0JacJNeyfB0Fbwf+7/f1clgv2SGtwVTEyMQkSJSJgn4BeSPZtnCk7asEEZobAgaZMCqhSdtnonSWHMvHb3cpjLmImJiZGJyLm1YoQkSUSC27LEosv2VpIofTxMwDtpOrPL7w2MDA2mnpSzrUJo1IeJiZGZNN5IMJwVJyJB9mzelrlzs25RgcO/t2xvpWlBgZmZW2nmU6tbUKy/pF0MNG2A0TvM6CtJISR7Nm+bISR7xrazZyy+csoiPn4e2balUL9MHrKMm5lBC0frB/tJskzM2blJPI3UmJjUxILJh2Z8epGBie1OBZhSSIL4RaQqQQnaUQeFBCUEf6Vc171SZ4d3kreV15al8+dMfYx6MTFpkF4SlKkZfwMikkVIgqR5U10vCUpuIryTJlv5vXRveth0L81iYlIDvfhgBgkTEg9/aCksrOUnw6tOS6FnBCWBuu6hOirgsvqA/OJhnkjzWAd8wzQ9niAN/oosb6ZWkDpF5ZFtW2rtmC+TXsuSatLWOj2RgTlzevaeqgoTk4rpda/EH97KIiRR/SELR5dlFpI82V1B6sj2Ojg+lj+zK8W7TqpueLTFK8nzgiyjeSzMVTMDE9unPh5tE5ywJIG0QpKUqZXU6V41vRj2Ct4rVSdwtO1+jMNmM24PJiY10gtpk6HTZsRM0likg9xIT9i900uVPtiI965jYtIS2lIx+D2nYHgrzRiOpBBW3r4SE63mSNsP0kt9O0b5mJgYU0wJmjeOJCAkQaIq+CjBqDuTK4peDHVBiNBTXyMkSSjKEhKb3LF3MTFpEU17J35PxD+GpIzKty1C0iWqCJvmuQfTCImFuKIRkbUi8hMR2SwiF8eUe7GIqIiscb+PisjDIrLB/VxZn9UzsWyuhoh6AJtKFfYqkYPjY6nEIynsVIV4lJHV1TUGJrazYHhZLffM5PCyWqae7ydEZBD4IPAsYCtwq4jcoKo/DpQbBv4E+H7gEPeo6upajE2gMc9ERI4TkW+JyF0islFE/tRdPyIiN4rI3e7fo5qysRcIm6YlT7aP55X0agio7whMaVNnuMv7pMW8kljOAjar6r2qegC4BjgvpNy7gb8B9tdpXBaaDHMdBP6Pqv4a8BTgdSJyCnAxcJOqngTc5H7vHHEhirQVQ1K5TMeJeT97F2mlaKad9DFQrqxwV9Nh1g6zRETW+z4X+batBO7zfd/qrptCRM4AjlPVfws59vEi8kMR+baIPK1809PTWJhLVceBcXd5QkTuwrmI5wFnu8WuBm4G3tqAiT1BXEWStvU4MLE9dXirSfoixJVi8GJU+SIhUhOSjMyem2WA6k5VXROxTULW6dRGkQHgCuBVIeXGgVWquktEngj8i4icqqoPpjWsTFrRAS8io8AZOPHAZa7QeIKzNGKfizyl37VzZ12mlkoR72RqZHrC8VN5LxbeahdxHkoLJ4Uskz4cfLgVOM73/Vjgft/3YeA04GYRGcOJ4twgImtU9RFV3QWgqrcB9wCPq8XqEBrvgBeRBcAXgDeq6oMiYUI9E1VdB6wDWH3mmZpQvLXkmXtpmkDEVTwpWrieV9J2qvBKWj1nV5b3nBT0TswraZRbgZNE5HhgG3A+8DJvo6ruAZZ430XkZuDNqrpeRI4GdqvqIRE5ATgJuLdO4/00KiYiMhtHSD6tql90V28XkeWqOi4iy4EdzVmYj7IezriKIWxk+gx2bWMAIrN9/F6JPzurbeGkttnTVryGSdT9F3kPGI2hqgdF5PXA14BB4CpV3Sgi7wLWq+oNMbs/HXiXiBwEDgGvUdXd1VsdTmNiIo4L8jHgLlV9n2/TDcAFwGXu3+sbMC83eR7OsBBFXOqw0U7Kfn1vJlowUWRZjAwN9tXgRVX9MvDlwLq3R5Q927f8BZzGeCtoss/kqcArgWf6Bt08D0dEniUid+PkXl/WoI2lkTWmHTfCOZVXUjFVewsLT1w59THSk+U+60rDpA/7WVpJk9lc3yU8kwHgnDptKYuoh7PM1M0mhWTuylXTOuoXnriy1Dmz6haOMvtLGvVKPFzvJHi/+b3cLN6Jv1wR4cnjEfWbd9IFGu+A7zpFhMQfA29DeGLW8tHSsr6a9jha2/FeAWGj1uPEIexe2zswr3ZBMXoLE5OS8D9oZaZp+gWlaa+ExStTeydNi0Uc/SQkYWQVEv+2OgXFvJPewsSkZLqS7x/EC+OEeSdtFo4gVQhJK0JcHgkd8XGp6GkqehMUIwoTEyMRrwKeHF7GADP7TnqBfvdGkshSwXtl6+rA9zrYw0Slsc73wdnZZiroA0xMjBkkzc5bZt9J1ZiIzCTPQNkw8o5bydt/Yllb7cbEpMtU0HLqFe/ERCSesgQlSNoMMOuQ7x6tmJvL6B32DsxjcnhZu/oJfMxduapWIWnldWh4DJLH3oF5sYLRlXEuhoOJiZEPN7OrLdQtIl2gDckied69Y7QTE5MeoA0PPaSbhbgJTETaTZpwVhvvKyMbJiY9QhumUGHXNgYmtrfDFpemhaQXZlyOoy0NFTBB6XVMTPqYNBXJ3JWrpvcLuK+LbcPLtJoWklaTUeyrrsits737WDaXkZo2tcJNSMqjqsyuPFiWV+9iYtIrNBRW8rySNngiHiYkBfHfS276+MDE9sj33pRF2tHzJii9iYlJL9DwfFxt6R8BE5LS8f63LRvN3XpBGZzdGm+uLVifSQl4ra1KOjOLVOS+fcNahEkPw6zlo/bAGIaRCvNMSqJNWTHTcF/dGyYKewfmcaQ7XmRhYNuUV9Iy2vTe9lYOWPRo4f8uC632SoxQTEwKMjU1fBWUGF6KiolHjWY3rySe1gpJGhEJhi5rEp60GWMmJL1JpJiIyInAMlX9z8D6pwH3q+o9VRtnlETMtOSedzLjRli8svB041XRtHfSSiFpuSeS5j4yEWkHInIacAow5K1T1U8m7RfnmfwD8Bch6x92t/1WRhs7xzSvpOyWXtmd3m64K8o7CXaemVcSTmeEpMbEiiQhMRFpDyLyDuBsHDH5MvBc4LtAopjEdcCPquodwZWquh4YzWNoFkRkrYj8REQ2i8jFVZ8vK9MekOBD6Q7sm/obc4KxAAAgAElEQVRkpebsKW/yRv+n7Q94E2nKrROSxSuLNVxq8GZMSJJJqutEZK6IfM7d/n0RGfVtu8Rd/xMReU4J5rwYOAf4har+PvAEYG6aHePEZChm2xHpbcuOiAwCH8RRxVOAl4rIKVHlD01WaU00qftKsohDQ2m43gyvoTO9tjSEUregtGnQZmn/k4r/t20Ti7a9tTFlXXch8CtVPRG4Arjc3fcU4HzgVGAt8CH3eEV4WFUngYMiciSwAzghzY5xYnKriPxhcKWIXAjclsvM9JwFbFbVe1X1AHANcF7cDnXeJLFeSRRJ5fJ6MRXjn3I+qp+iyf6LtgykrJWWinub2b3/UOuExCVNXXcecLW7fB1wjoiIu/4aVX1EVX8GbHaPV4T1IrII+Cecev4HwC1pdozrM3kj8CUReTmHxWMNMAd4YX5bU7ESuM/3fSvwZH8BEbkIuAhgxbHHAYcFpY43suWa7LCFYpGaxSthfKxpK0Kps0P+4PhY+8JdLScukaOOwYktEJElIrLe932dqq5zlxPrOn8ZVT0oInuAxe767wX2LdTSUNXXuotXishXgSPDujvCiBOTjwMvx+kfOc1d9++q+s28hmZAQtbptC/OP2MdwONXnzltW5WiUmkqcEvZOzCPBb4UYr834J8IskkvoekML6OdVCUkh2QgiwjuVNU1EdsS67qYMmn2zYSI3KSq5wCo6lhwXRxxYa5PAF8Dfh24UlU/UJOQgKOwx/m+HwvcX9O509HDXob3QqIsab97B+bB4pUzWuWzlo+Grm+CusSs0b6THg1xNfHGxRZ4JGlIU9dNlRGRWcBCYHfKfVMhIkMiMoLjRR0lIiPuZxRYkeYYkWKiqtcCZwBH4sTR3iwib/I+eQzOwK3ASSJyvIjMwelkuqHicxox+B94zwMIegJt8AzqFJRGRKWHGzFpX5LVxrFNFZKmrrsBuMBdfjHwTVVVd/35brbX8cBJpOzfCOGPcLozTsbpJ7nN/VyPkyCQSNII+EeBfTipYcNALXlTblzw9Tie0SBwlapurOPcnSZkapWiMeuD42NTN1Ee76SKyr/uPhSoOW04ZhBqL9K2jK86iarrRORdwHpVvQH4GPDPIrIZxyM53913o4hcC/wYOAi8TlVzuWOq+n7g/SLyBlX9QJ5jxI2AXwu8D0f9zlTVWpsLqvplnEEzRsn4+3wmh5clCspUS9HXKo6qrPtRUKABUSlZUOoapJrUGe9f7heRCavrVPXtvuX9wEsi9r0UuLSoDSLyFlX9G1X9gIi8RFU/79v2XlUNG8A+jbg+k7cBL1HVi+sWEqNCAmGStMkEceVyhXvcAXdxKcdFeWTblkbGotQW/urhkJeHJxhVhLbqyOrsEOf7li8JbFub5gCRnomqPi2PRUYPEDL1S1RLMHP2WtbWcg0px56gdNpT6SHsJVmtRCKWw76HYrMGG1OvbY19wF0Bimp1e+Mv8oRL6kotbiJ9uHJRKTnc1UQF3mcd7m1FI5bDvofSaTHZvf+QubpJJFRGebySpP6X4DHrzIpqajxK2z2Vut8Dbx5H63iCiDyI44Uc4S7jfo+bWmuKTouJkR6vMgkVj5Lfq1LVsXuBXhhB36bwUpts8XNosmfGsaRCVQu3uu21vcZUhZ7khRT1IKZlhbmfJsZqND2fV6smjPTR1pkdLAzWG5hn0kb8Yae6Wu5R4a6yvZKGBKRttN1DaatHYLQXE5MSKDUeHqzQa3yJUdx5yvJKDo6PNe4ZQDvm8mq7oLQJE7f2Y2GuttMjfQqpKsUe+S11UoqXZtfVaAHmmZRAqa3LFlYMWSq8Iq1Hz1Ooy3Npg3fS71h/SHcwMTFyMUNAF69MTi2NGaDor9TnrlzVilBYXVi4y+gCFuYyYon1Srx3kKcYXxJHmHdQl8fQGdFqoUebhHkl3cI8kzLwOs178IGOI0pIvHeY+Mky4M1CS0ZdjAwNdmo8SJsxzyQvXmXq/p0cXja13IXU16TfMDm8bNqnikybfvJO2tARH2wQVOk55Dm2eTLtxjyTDETdzMGHsAvx71nLR0MrOL9XkkdAoo5rlETH3nVSBlV4J49OTrJj34FSj9nrmGdShD57aD0hyeOJ7B2YV+vcT71IaSKb00MJGwHfNm+gbfYYhzExSUkZN7H3rovGXvmaEb+HVZq3lfF98f0U6oKSBaVjfXgeeZ5Fm/C1ekxMUlBESOLEoxdEZdby0cOVf06vxMO8k3SUek+UIChd8QZMUKrFxCQHaSvEtJVCL4hKqWT0ToyCdNBD6YrAdYlGxERE/lZENonIHSLyJRFZ5Nt2iYhsFpGfiMhzmrDPT9hNW9Xsqq0WlIJeiceUd5Khvykq1DV35apSw2BtCXVB++6FNlbeFu5qF015JjcCp6nq6cBPcd85LCKn4LyL+FSc9w5/SEQK/fd7Lce8bZUIUJqQeHjHydp3EhwlH7UtL50e/9Iy76Sse2nB5EOtFLp+pBExUdWvq+pB9+v3gGPd5fOAa1T1EVX9GbAZOCvveYqm7qXxSqroA+ilTvq8eN5J1nBXnHDkERVvnzYKSRP//zivu45KO6/X31VREZEREblRRO52/x4VU/ZIEdkmIv/oW3ezG+XZ4H6WVmVrG8aZ/AHwOXd5JY64eGx1181ARC4CLgJYcexxkQd3BGVOz7q3jc/bVLJX4rF3YB5HlnrEw7RRGFpBy8ag7B2YF9tgy9NQ6+BU9RcDN6nqZSJysfv9rRFl3w18O2T9y1V1fVUGelTmmYjIN0TkzpDPeb4ybwMOAp/2VoUcKvRl9qq6TlXXqOqakcVLEu3JGu4q2ldSZquyqx5KXu/EaI46W//mpQBOtOZqd/lq4HfCConIE4FlwNdrsmsGlXkmqnpu3HYRuQB4AXCOqnqCsRXwuxnHAvfntWHp/Dns2HeAHfsOsHT+nLyHAdr7StNKqcgr8dg7MI8Fw8sYwLkRuyqaRWjcM62RsGdsYGJ77lBylV7KgUPKlj370xZfIiJ+z2Cdqq5Lue8yVR0HUNXxsDCViAwAfw+8Ejgn5BgfF5FDwBeA9/jq21JpJMwlImtxXLXfVFV/E+IG4DMi8j5gBXAScEvS8QYtwTkbaV7P26JwiNHfFBGUlrBTVddEbRSRbwDHhGx6W8rjvxb4sqreJzIjuPNyVd0mIsM4YvJK4JMpj5uJpvpM/hGYC9zo/vjvqeprVHWjiFwL/Bgn/PU6VW1FOtbk8LJGvJNKXgUcVaZlGT+GQ6neScv6TdKSV1B6oQ8lLoojIttFZLnrlSwHdoQU+5/A00TktcACYI6I7FXVi1V1m3uOCRH5DE5CU3fERFVPjNl2KXBp1mOmmcxt9/5DPdMRX7jyyFNhBASl7tagTQLZX/j7NdI01LosKDHcAFwAXOb+vT5YQFVf7i2LyKuANap6sYjMAhap6k4RmY3TrfCNqgy1AFEG6qpccwtJ4GVVhY5htAoT2b7lMuBZInI38Cz3OyKyRkQ+mrDvXOBrInIHsAHYBvxTVYa2ITW4NNryIpwiD34uIbHK3yiJuhpMWb0Sf9l+8k5UdRchnepuqu+rQ9Z/AviEu7wPeGK1Fh6mU2LSBNMehIJ9DiYkRj9QNG233wSlV7AwVwbStKBqDUd0SUjc39IvabBZ6ep1yZPU0uOZXZ2lc55JY6GuXdt6SkiCD2RfjqPpAUoXkYYbIE0PJmz6/F2mc2ISxBu4WCpNp9AWqBCiWnXB1Gdr/TVHU15I3f9za8B0i86LSVlM3fgRQhL0SvzTmaeZK6rqCiRNRdHUWJqmzts2uhrKKpO2NHIOHJxkbLd5OX5MTLIQEJKwsFbl78TI4ZVkeQCLPqx1dnKWMS6lSAVeVlizyyLihZyD47us8dA9OikmbUkRbgNVteTiRKMOQfEqYH9FnLVyL1qJFzl3GefPTI8nbLTFKzHC6aSYBCk6yWMYNohsJnV0bsZVwN62NP+bsivy1nsXDQlJ1IwTaUObJiC9g6UGR2D56IeJEok48cglLDHpwbOWj6ausJPKtb7iL5MMMxqkqdyLPBfBeyJJKExIeovOeiZRoa6o8Je/BZU0h5d5JfWRt+LP4qV0EgtpGTXTWTHJSlBkdu8/BEPhb4JrlB6a9bVI30kmEfGuR0imXd95IRVShbduKendoa/FxBt/EtenEvVq0SjmrlxVfUZXw/ivhycYjYhusPKMEZWepeaGQxPvDjEB6QZ9LSZVkfUd5Lla4TWT18tIyvqaQdr3qsRdh15+N0vDXmdUxV6GV9JYo8N3fqM6TExyUEY8vo5JHato8YVVBrVWEGmvQZOC0iNhyLSUWQmXISgmCu2ks2KSZpxJ5qlWAhVUHlGpa2bguoQk7X7BCqBItlfm8mWJSsdEIkjYPVO04g5LZskiKG0Vjv2PHmLT+ETTZrSKzopJWtKOQYnLiy+1k7eECquTMeiE6+L/zTP+T0W8lI4LSBxlCIn3N0xQjG7R92KSmTJbuxVUVG0VEr93ktkrySAk/u/TRCXr/63PRKTs+yYsO7JXXplt5KPRQYsi8mYRURFZ4n4XEfl/IrJZRO4QkTObtC+WDr8ad+/AvGZbjv5rU/YMyUnXvuX/myykFYi6GiC79x+yaY46TGOeiYgch/NOY38e7XOBk9zPk4EPu39jOTSZro8kjB37DkSGurzWlBfjnRFKKanSKePdIkUqhCjhyCIoabyNzBlhJYb8QkNfPUae/3GR+yJvBl/Ss+jfnsdbMUFqJ02Gua4A3gJc71t3HvBJVVXgeyKySESWq+p4FQZ4ne9xguIRDNFU2ZpLOnYZYwHK9jzShrBSh7hKruzbOM19W0OSdWLC0B0aERMR+W1gm6reLiL+TSuB+3zft7rrZoiJiFwEXASw4tjjMtuQNosrOIV2WCWcJzMprjJPOl4VXkhZ+I/fttkDIr2UCs/V62T1TsLEIc3gYCMcERkBPgeMAmPA76rqr0LKXQ483/36blX9nLv+eOAaYAT4AfBKVS35bYEOlYmJiHwDOCZk09uAvwCeHbZbyDoNO76qrgPWATx+9ZmhZbKQdMMHHxK/e15VK9/D8vLLpaiX0hWhiMO7PpPDywpNi+NvtAUbcCYuqbgYuElVLxORi93vb/UXEJHnA2cCq4G5wLdF5Cuq+iBwOXCFql4jIlcCF+J0H5ROZWKiqueGrReRxwPHA55XcizwAxE5C8cT8bsZxwL3l2lXkkeS9oav0j2PSqPMIiptEJCmRzzHESYIYQLTD8IRpCzPLeuz5idOaEp/DXe7OQ84212+GriZgJgApwDfVtWDwEERuR1YKyKfB54JvMy3/zvpNTGJQlV/BCz1vovIGLBGVXeKyA3A60XkGpyO9z1p+ksenZys7AbLctyyWlpRHZRtEIgyKFNk0lT2aT29rghHpulrQsjrufnv26LPY8cEY4mIrPd9X+dGVtKwzKsDVXVcRJaGlLkdeIeIvA+YBzwD+DGwGHjAFRk43G1QCW0bZ/Jl4HnAZuAh4PebNScbZT4AnjAVzXxpmqBweBVdcH2eCixv5d82j6muRkKW/iz/tU1jX5yQbNmzP3bfVQuHEo/fNh559BBj21OPgN+pqmuiNiZ0CSSiql8XkScB/wX8Evhv4CAZug3KoHExUdVR37ICr2vOmmSSHgxI93BEHcfb1/9A5hWWPKG4KgQrLvU4b6WeVkiS0p6rFpW2epNpvbUiQpLmWUlTLup5Snv8thPVJQAgItu9jFYRWQ7siDjGpcCl7j6fAe4GdgKLRGSW652U3m3gp3ExKYMDhzT2xirS8slzwxa5yf37BoXFH0arqs8mzcvDqiIu26qqEFRbK/u6CROXrELip8yKviuikZMbgAuAy9y/1wcLiMggsEhVd4nI6cDpwNdVVUXkW8CLcTK6Qvcvi06ISRJV3Ixju6tt0Y6OzJuyO05U6iJOvIoITZh34oW82jaWpp/Ic+28+zP4vGV9VkZH7P/m4zLgWhG5EGeA90sARGQN8BpVfTUwG/gPN6HpQeAVvn6StwLXiMh7gB8CH6vK0E6IyYGDk7WdK+zB8M8eevLy4VLPEycqEC0sUf03VQhRXNp0XrrSGd5lwsJbfiHJ2+CK288vNFU36NqAqu4CzglZvx54tbu8HyejK2z/e4GzqrTRoxNiAtE3VlmtnCQRiVuXB0+UxnY/NPUbgqIC2Tv986ZjZiFr304VHeL+lnXWkGBZIb22je7O8ru8/0echzIyNDj1G73XOaxaOBQbCYh6PtI2wvpBQHqVzohJFGlbOVn2q+M9BpvGJ0IFBcL7VcqgCqFJKyxlCkpZU6d3jbT/i6qSEuKem7BtZXn5Rj10QkyCL6qpspUTdtMnpQiOLsv3UAQFBWYKYNH+oLRiVIbQJIXDigpKmIh0VRiKEjUlvHf9vSSIBcPLMnsnYeRpgNnLp3qLTohJkLJbOVE3ddo882C5OHHxlx1dNjxNUGCmABYN44WJUVZvJ2//THDeM8iWspunkz5rWLDsPqYyxiKVGY70X/tpQrJrGyxeycDE9kyCkpaoZydvw8tonk6IySOPJt/IWWO1ca2iDIOVEvf3PzzB445tn5gSFI+gvXljyHEiVIbAQLpEAcj3alev5ewJSprQVp6KvI0jscsMRwav/bS0bFdQshLsN/Hfv0nPTpaGVxJFn1MjG50QE8jf0snqSoed54Ff7kvcb9HR81MfL7jd/xvK8rqyejhFBSZJWOJe7RonKnGeSZnTe/QKYb8zi8BUOVV/nsq9rYJw6OBkque+n+iMmEQR5QUUOY6ftDfUA7/cFykoac4dZ3tWQQwTnzAPJ4/A+IkSm6jJNNPOSRZ8UZm3PW768yx2J1FW0kNaO4qeL2ny0uBL4IApjyTM6wsKf/C6p/1dUc9O3ufEaJZOiIm/lRB3I2bxXpJaRHlaJUUFJY4sQpnWuymabp00ZYxH0gj/OI8lrZD08ojsssKOHlHXe0pQfINGw4TE/zcvcc+Pf1sZwmIeRD10Qkz8BG+cNDdjVle6yM1ZRFDiSPsbokQnKDBxobO86dYewcoxaTBmkscSrNjSTjRYdjJDnZQxfVBQVIKCEiSPgHj3lf/+zPL8mBD0Dp0TkyB5xCXtsYocJ60dSedcdPT8VF6ZR5johAlM3r6ZNAkBUenNUYMx44Qly2y1SbZVNSAuSaTyhBjjSOsRevhfW+0XlCoIu58f2DH9Xlu01DK6epFOiMmhg5PTbsi4mzGucg5WxmnEI/ggzDhmhC1pBCXV+X1l0ghPGGkzaIqOXp46n6/yjBqMCclzkqWdrbbpUdN5zl+2wEC4aHsEBaUofvvjvOaw5yfts5yVpGfVKEYnxCRI3psxk/ud8sZ8YMdEIUEpk7ReWtb0zCwJAEmpzXEj/dNkZKWZG6roYLimRmaXNWVQ1jBjFFHJDZ6dwescvP/SPENhZdI80yYc9dNJMfFTRSsn642aR1CShC2NDUm/N6+4QLFR/UHSDMqMm0ImjSdS5mjqIscKE6KiE4XmCS36iXvtQRJR137T+ERkP0nRit6Eop10Xkz8FI3NFrmJkwQFDlfmZQhJWLmyxAXy5f9nCZ8Fp5GBcGGZZlOgUm3jdBxJNiVtz+sVpfVq8mSq+Y8d7HBPIyIT92+e9n14xYmZbTCapxNicujA/mk3ZNqbse4WTpygQDl9NFn2LVNc0pAlvTksuyxt30PR6W/KooqpQeLEpkyvJigySdc+bJS7d/+kFZHg+rJFJep8Rjl0QkyC5BGWukgSlKR9y7bFT1ZxCVKF2HgVclzqcpVT3xQh67mLik9aT6xIVl6aax3njWSp0MPKZn2eqxKQYNKP0aCYiMgbgNfjvPj+31X1Le76S4ALgUPAn6jq14qcp40udB5Bydq6C5Lmd+ft7Jzav6TUab8oRSUDlCEgYfbmzbArI5Eijd1leDtlhP+SphXy30tlVujmXbSXRsRERJ4BnAecrqqPiMhSd/0pwPnAqcAK4Bsi8jhVLW0e8Spd6CzHLOKh5CHv764qTTP2nDEjoPN6GlmmvSnz+GVn61X1uoMi54zL0rLKv39oyjP5Y+AyVX0EQFV3uOvPA65x1/9MRDbjvHLyv8s2oKiohD0kWcNraQWlTHe6SAgwyo4qRSbv1BptGTmdx/spQh1hvahrWySkZfQ+TYnJ44CnicilwH7gzap6K7AS+J6v3FZ33QxE5CLgIoCBoYW5DUkTl83zUKQVK+8BjMz0qjAuW1bfUtk2JmW99TpZf0eTEx/mzSw0Iek/KhMTEfkGcEzIpre55z0KeArwJOBaETkBkJDyGnZ8VV0HrAOYtXBlaJm8VBHjTSsqTdGmvqWifTdpjlclZXtqacUnj+hkFbaka2ki0r9UJiaqem7UNhH5Y+CLqqrALSIyCSzB8USO8xU9Frg/6VyTjz7CxPg9AAwvf2wRsysja59K05SRSVMmaaataVqQPcoWw9TnrcBzy3JNM2Vq1fi8eufqRUTkJcA7gV8DzlLV9SFlhoDvAHNx6vTrVPUd7rZPAL8J7HGLv0pVN1Rha1Nhrn8BngncLCKPA+YAO4EbgM+IyPtwOuBPAm7JcuCJ8XtMUCqiDO+lqrTttghJFE30N+Ulrffx+Gf9xvT9lg5z34YfziwfeCb9lXtVotLLAhLgTuB/AR+JKfMI8ExV3Ssis4HvishXVNXrMvhzVb2uakObEpOrgKtE5E7gAHCB66VsFJFrgR/jpAy/rsxMrjZQtqAMrzixsdBCmhBenG1tCq0lXcOqbMs7UWidNvjx/89nZI4tG+aBHTPvx6BQDC9/bKWVfYeEBFW9C0AkrAdgqowCe92vs91PqaH/NDQiJqp6AHhFxLZLgUvrtcgoQlliVkVorc22paFNHpf/GoyeegxrT53+zpOx3Q+xYekwE4mB6cMCU3bF31IhWSIi/vDUOrfPtzREZBC4DTgR+KCqft+3+VIReTtwE3Cxl0VbNp0cAd92qvBOvON2jTb/pjbPtFA2/t963OozWLt6BWtWHDmtzKqFQ2w67ZhQ7ySKssJbdYtIcAqnBHaq6pqojXHJSqp6fZoTuBGc1SKyCPiSiJymqncClwC/wOlKWAe8FXhXWsOz0EkxabLfJG3FXsXgyS6LSttpU8iuLKLuo0VHz2d0ZN6M1ymPDB3B2lOXMbZ9gh/lvAfDRCHqWW6pF5KZuGSlHMd6QERuBtYCd6rquLvpERH5OPDmss4VpBNiMjB7bu3iMbzixNh49qKlZ9RoTXvOnYewTttep9fEJW0DJMorKYssz3HeZ353JblMzSEiRwOPukJyBHAucLm7bbmqjovT6fI7OB36ldAJMWmCRUuHWX1amGdqQMaR2KvP6KSg+GlrSCyLFzu84kRWn3YMTz9+ZIZXArBp18N8deN2xjb+okwT+xoReSHwAeBo4N9FZIOqPkdEVgAfVdXnAcuBq91+kwHgWlX9N/cQn3bFRoANwGuqstXEJAfeQ9XUG/d6gSzXZtP4BPd1rLUYRxvG8OQJhS5aOszJy4cj38K4Zc9+xrZPWJi1RFT1S8CXQtbfDzzPXb4DCA1HqOozKzXQh4lJDrxMlrB3aRv52NAH3kkceSrgOqdjN6/ESMLEJCNefv2qhUOp3pNtJPP040cyZwEZ9SZaeA2okxcfEbr9Oz/bzYY7f2H/vz6mE2IyOGeotjBB0kNl5GEOJy8fZsPSYZw0+XiswppO2L1ftqfjNaDC2LTrYTaNT7RqTEzcb9ldox39RCfEZM7QLEZPraczvMpMln5lZGhwyjvZkCrXxDwYP4uWDk+b5NGZoyv7NYp7huIaUG3sKwleEz8/r9mWfqETYjJ39mAlLwUKIyy/3ijOyYuPmDGiOoqx7fNzj2PoGsetPmNGVuEYzsj5LFPthE6P4qOX+ge9/p0ovl3COfyTyxoOnRCTI4+YnboiKop5JdWxauFQqv/jV4GxBuckaxOLjp4/c1qT5cNc88t9mQQlbHoUP70U1k36LVfXaEs/0QkxmT97sLZK3ryS6khbYRUdZd0VvAGEQa/Bm9bk5m+m68PwH6fXk0r8CTJGvXRCTAYHrJLvCmkEZen8OaxdvYKxjf3tnXjTmgSv2e79h6YSGtJ4J1HH6UUsQaY5OiEmRn8xMjTI6Mg8Rk89hrEM+7Xx3SFFiEoGCUtoiJtip2hSyaqFQ4wuG2asxIzKvP8r80qaw8TE6EnWrDiStatX8NWmDWmI0WXDkQMI4XBCQ5ppbeKOk4Ys56oa80qaw8TE6Ek872Tt6hVNm9IYSf0bqxYOpbo+ZfSTpD1X1ViCTHOYmBg9S56Ko0shkCRvwmmhjxQ+ThpOXnwEW/bMK3wcj7z/J+s7bQ4TE6NnGRmqL4uvbaStNOvMzurX/4XhYGJi9DTWEo2nzutj/4v+Rpx30fc2IvJL2jVLwhJgZ9NGhGB2paeNNkE77WqjTRBt12NU9egiBxaRr7rHT8NOVV1b5Hy9QCfEpG2IyPq4dz43hdmVnjbaBO20q402QXvt6ioDTRtgGIZh9D4mJoZhGEZhTEyqYV3TBkRgdqWnjTZBO+1qo03QXrs6ifWZGIZhGIUxz8QwDMMojImJYRiGURgTkxIQkTER+ZGIbBCR9e66ERG5UUTudv8eVYMdV4nIDhG507cu1A5x+H8isllE7hCRM2u06Z0iss29XhtE5Hm+bZe4Nv1ERJ5TkU3Hici3ROQuEdkoIn/qrm/6WkXZ1fT1GhKRW0Tkdteuv3LXHy8i33ev1+dEZI67fq77fbO7fbRGmz4hIj/zXavV7vpa/od9jarap+AH502pSwLr/ga42F2+GLi8BjueDpwJ3JlkB/A84CuAAE8Bvl+jTe8E3hxS9hTgdmAucDxwDzBYgU3LgTPd5WHgp+65m75WUXY1fb0EWOAuzwa+716Ha4Hz3fVXAn/sLr8WuNJdPh/4XI02fQJ4cUj5Wv6H/fwxz6Q6zuPwG0KvBmZQBYYAAAQJSURBVH6n6hOq6neA3SntOA/4pDp8D1gkIstrsimK84BrVPURVf0ZsBk4qwKbxlX1B+7yBHAXsJLmr1WUXVHUdb1UVfe6X2e7HwWeCVznrg9eL+86XgecIyJSk01R1PI/7GdMTMpBga+LyG0icpG7bpmqjoNTSQBLG7Ityo6VwH2+cluJr7jK5vVuuOEqXwiwdpvcEMwZOC3b1lyrgF3Q8PUSkUER2QDsAG7E8YIeUNWDIeeessvdvgdYXLVNqupdq0vda3WFiMwN2hRir1ECJibl8FRVPRN4LvA6EXl60walIKylWFee+IeBxwKrgXHg75uwSUQWAF8A3qiqD8YVDVlXp12NXy9VPaSqq4FjcbyfX4s5dy12BW0SkdOAS4CTgSfhzL//1jpt6mdMTEpAVe93/+4AvoTzsG333Gj3746GzIuyYytwnK/cscD9dRikqtvdimAS+CcOh2Zqs0lEZuNU2J9W1S+6qxu/VmF2teF6eajqA8DNOP0Oi0TEm3ncf+4pu9ztC0kf6ixi01o3VKiq+gjwcRq8Vv2GiUlBRGS+iAx7y8CzgTuBG4AL3GIXANc3Y2GkHTcAv+dmuTwF2OOFeKomEKt+Ic718mw6380GOh44CbilgvML8DHgLlV9n29To9cqyq4WXK+jRWSRu3wEcC5Of863gBe7xYLXy7uOLwa+qaqlegERNm3yNQYEpw/Hf60aud/7hqYzAHr9A5yAk1FzO7AReJu7fjFwE3C3+3ekBls+ixMGeRSnJXZhlB04bv8HcWLfPwLW1GjTP7vnvAPnIV/uK/8216afAM+tyKbfwAlx3AFscD/Pa8G1irKr6et1OvBD9/x3Am/33fu34HT8fx6Y664fcr9vdrefUKNN33Sv1Z3Apzic8VXL/7CfPzadimEYhlEYC3MZhmEYhTExMQzDMApjYmIYhmEUxsTEMAzDKIyJiWEYhlEYExOj53Fn2/2ZiIy4349yvz9GRB4nIl92Z4u9S0SuFZFlTdtsGF3DUoONTiAibwFOVNWLROQjODM5X4EzpuBNqvqvbrlnAL9U1TsjD2YYRmZMTIxO4E5DchtwFfCHOJMkvgI4W1V/r0nbDKMfmJVcxDDaj6o+KiJ/DnwVeLaqHnAn/rutYdMMoy+wPhOjSzwXZ+qW05o2xDD6DRMToxO4r2d9Fs5stn/mTvi3EXhio4YZRp9gYmL0PO4MsR/Gef/HFuBvgb8DPgP8uog831d2rYg8vhlLDaO7mJgYXeAPgS2qeqP7/UM4L0g6C3gB8AYRuVtEfgy8iubeLWMYncWyuQzDMIzCmGdiGIZhFMbExDAMwyiMiYlhGIZRGBMTwzAMozAmJoZhGEZhTEwMwzCMwpiYGIZhGIX5/9toWN6WfSzMAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ds.Eta[0].plot.contourf(levels=30)\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXu8HmV1778rCUkM2STs3AwJuIFAKSACRmq9IiBG9IgetQdBCx4sx6qt1mM1VI/lIFToBeqpF8wRJFYxINWStqh4A+tBhaAECaAE3EJCDCSBsAGTsJN1/phndmbPnpl35p3rO+/6fj7z2TPzPDOz3tkzz+9Z67mMqCqGYRiGATCpbgMMwzCM5mCiYBiGYYxhomAYhmGMYaJgGIZhjGGiYBiGYYxhomAYhmGM0QhREJF1InJi3Xa0FRE5S0RuqtuOshCRC0Tky3XbYRhtoBGioKpHqerN0KwXXEQ+ICIPisiTIvKIiFwuIlMC6S8RkdtEZERE7hKRl4WOnyci14jIEyLyuIh8JZA2TUSucuf+rYh8MHTsySJyn4g8IyI/EJHnpbR5SEQ0aKeqfkVVT+3+TpSPiBwrIne433uHiBxbt01xiMhUEbleRIbdvT4xlP4q9z/bLiLDKc435PI/4/7np4TS/8I9I9vdMzOtinOFjlsnIk+Flp0isifN8d0iIoeLyA0i8piIbBORb4vI7wXSzxGR3SG7Tkw4X1fvVT/RCFEokmBhWAD/BhyvqvsBRwMvAP7cXWcQWA38HTAb+Fvg30Rk/8DxXwd+CzwPmA/8fSDtAuAwl/Yq4MMissyde6479n8Bg8Aa4NoCf1ejEJGpwA3Al4H9gZXADW5/U/kR8Ha8/2+Yp4GrgL9Mea6vAj8H5gAfBa4XkXkAIvIaYDlwMjAEHAL874rONYaruM30F+C5wIPAJ1L+xm6Zjfee/R6wALgN71kJ8uOgbX4FM0y/vVddo6q1L8AwcAqwDNgFPAs8Bax16bOAK4FNwEbgImCySzsH+H/A5cA24KKSbJwDfBf4rNt+PbAulOdXwLlu/VT3uybHnG8jcGpg+xPAKrd+HnBrIG1f4HfAESnsfAhQd/+eAv7Q3aMfBfIo8B7gfmDEXftQ4MfAk8B1wNRA/tcDdwJPALcCxxR8b09190NCv2NZTP6DgVuc7d8BPg18OZD+NbzCejvwQ+Aot/9FwGZgSiDvm4E73foJeAXFky7fZSls3wCcGJN2CjDc4fjDgZ3AQGDffwLvduvXAH8TSDsZ+G3Z50rxu1cBNwGTynjfEq476J7fOW573LPd4diu36t+WhrlKajqt4C/Aa5VT/Ff4JJWAqPAEuA4vELkXYFD/wCv1jIfuDh8XhE504Vw4paD4mxyxz4JbMHzFD7vJ7llXHY8jwLgxcAvgZUislVEbheRV7pz7g8cAKwNHLsWOMqtHxVMU9WngQcC6Um8wv2d7e7hj2PyLQNe6Oz8MLACOAs40P2Gtzlbj8er9f4PPGH8PLA6Luzgwmhx9/mzMbYcBdyl7k113JXwe68B7gDm4gna2aH0b+J5YfOBnwFfAVDV24GtwKsDed8O/LNb/xTwKfU8w0PxxLFwROTfRWS52zwKeFBVRwJZYp8Ft75AROYUfa4M9v858FLgTFWNDB+JyMs6vHMvizouBa/AE7KtgX3HicgWEfmViPyvhGhBnveqbygy1FIKIrIAeC1eIfc74GkRuRxP9f0C+hFV/Se3Pho+h6peg1eQZMY/VkQOA/4YrwYJXo35ABF5G3A9cCZeQTLDpS9mr3i9E69GeoOILAGe4/JsD1xqOzDg1mcCj4VMCaYXwaWq+iSwTkTuBm5S1QcBROSbeOK7EvgT4POq+lN33EoR+Ss8MbklfFJVPaYLW2Yy/l5AzO91Av4i4BRV3Qn8UET+LWTDVYH8FwCPi8gsVd3uftPbgW+6EOBr8Lwm8DzUJSIyV1W3AD/p4rd0RFVfH9iM++2LYtL99QFga5HnSmO7iLwYr+J2irtHkajqj/BCP4UhIouBzwDB9rcf4lVifoNXuF+LVwZ8MuIUVbxXPU+jPIUYngfsA2zyaxl4YjA/kOfhso1Q1fuBdcBn3fZW4HS8B3QzXs37u3jhBPDc0mFVvVJVn1XVVc7Ol+KFdQD2C1xiP7xwCC49mBZOL4LNgfXfRWzPdOvPA/5nsJaH500cUKAtWX7vAcDjrpbn8xt/RUQmi8glIvKA8/CGXdJc9/fLwH8RkZnAHwH/qaqbXNq5eCGY+5xnFyxwy6LTbw+n++tR96bIc03AxeS/BpyvqqUIZqjB+KDA/nl44arPqupX/f2q+qCq/lpV96jqL4ALgbfEnL6K96rnaaIohKdtfRgvTjpXVWe7ZT9VPSrhmHGI1yUz3HMi8uHrwBQ8b8C7qOotqvoiVR0E3oHXGHabS74rzi5VfRyvfeQFgd0vwBMd3N+xNBHZ1113HZ0petrbh4GLA/d+tqrOCL6YQWJ6qfjLFTHXWAccIyLBcNwxRP/eTcD+7p74BP9/Z+KJ9Sl4bVFDvmkAqroRr+3kTXj/Mz90hKrer6pvw6twXIrXSBu8ThmsAw4RkWBtNfZZcOubQ+GTMs41DhGZhOdt/7+AV56U/+Ud3rmXRx2n4xuMH3Ln2h9PEFar6oTwcPgUTAzr+uR5r/qHuhs1XBh5GM8dBXg3Xs+OSYH0G/DivfvhCdmhwCs1Y0NTF3a9C5jv1o/Ee3guC6Qfh+fF7Af8I94L46cNAo/jxbsn49VetuGJG8AleOGX/YEj8Aq7ZS5tHp5b+2ZgOl4B9ZPAuS8Abo6xeQawGzg8sG/cPcJ7cZYEtn8EnBPYvgj4gltfiicMf4D3su0LvI5AY2YB93kqXm3//cA04H1ue2pM/p/g9eSaCrwMr2H4yy7tPXiN4vs5Wz8b8XvPAn7hjts3sP/twDy3fgqwA5geY8M097/ZgBcmnI5rKHfP6HS8sOdv3Hrkbwn9nul4YvVEwI5leI3mR7pn5fvAJWWcC7gZuCDmvBcC9wIzyywLIq67H15F69Mx6a8FFrj1I4C7gb+OyZv4Xtni7lPdBrh/1jB7RWGOK6QeB37m9s0CPudewO14Xe7OcGnnUJ4ofBEvrPK0s/HvgoUEXve/7W65FicggfSXu8LnKbxeLS8PpE3Da8D1e7p8MHTsKcB9eKGcm4GhQNqVeLX3OLsvxIudPoEX+x93j8ggCm57GXC7O98mvBBCYaLgrnEcXuPx7/Aah49LyHsIXq+apwj1PsILe92AFxL4DV47UPj3znD3fWXovF8GHnXnXQe8scMzq6FlyKWdGJF2c+DYbwJ/Fdgecv/j3+F1TjgldC0/RPmkeyanlXSuB4BXx/zePXge+1MRy0FlvH/uume7+/d01DXxBNB/Rx/Ee/b3CRy/Djgr5Xt1BXBF1LF43mipv7Upi1+zMXoIEbkTOFlTuP1GNCLyAPA/VPW7ddvSBFwj7tdU9Q/rtsWoFxMFo+8QkTfjhQ4O15gulYbRrzS+S6phFImI3IwXU3+HCYJhTMQ8BcMwDGOMJnZJNQzDMGqiFeGjfWcN6v7PXdQ5Y4DnDoyfpWHb754t0qRYBp+zz4R9cdeOymsYxl7i3p1do/GRwY2/unuLqs7Lc90D5Tm6g3TRxy3s+raqLstzvSpphSjs/9xFvO+Kr2c+bvmrDhtbX7V2Y5EmRXLGC6KFK+7acfkNw/CIe3eGtz0Te8z5Jx3+m9jElOxgD29mYaq8n+c3czvnag4WPqIaQUjCCn/DyE43gmB0phWeQrdc8oP7x3kLSXR60IYGZySmW8FvGEYv0PeeQlJtI7h0Im/txETDMNJjXkJ59LWnEEWeh2p42zORHoMV+EYvc+LlE2ZIj+Tmv3hlyZbs5aBZ03lo+45x+0wQiqHvRaFJD5KJh9HLnHj5LZUIw63DE2d3adJ73Ov0vSiUjRX0hlEcUWEjE4Ri6StR6NQYDPkesDTnN4w2U5W34GOCUDx90dA8NDgjdYGdNq+fz89rbQmG4ZG2DSIrQS8hbQcQIzut9RTy1tqDx/sPX9w5zUMwjOowMSiXVopC0YV0N+czL8FoK0MLBhjeXO1njX0voSmCMG2ScOhzpqbL/HTnLE2ilaIQhRXShtEd4TaCc675WS3CEOS+TXuvfcTCgYScRlZa16ZgoRzDqIahBdUUxpf84P5KrmN4tEoU4gTBvATDKI6rzzx+bL0qYYjj6jOPTz1VjZGO1oaPDpo1vW4TDKMvKDOUZF5C9bTGU4jzEl4yNKdiSwyj/QS9hSbYYN5CcbTSUzAvwTDScc41P5uwr4gCP8sgtrA30G+NyCJyFfB64FFVPToi/XTgE8AeYBT4gKr+yKXtBn7hsj6kqm/Ia0/rRCEoCOYlGEY8UYJQxbFBkgTB3z5i4cC4rqj+viT8yEFTurB24Grg08CXYtK/B6xWVRWRY4DrgCNc2u9U9dgijWlF+GjqlFb8DMNoBFkL/Lj2hCwjm+/bNDJBEDoRFJRebntQ1R8C2xLSn1JVdZv7AhqXtwha5SmYl2AY6WiilxBH0DMIC0encwwNzmiCtzBXRNYEtleo6oosJxCRNwGfBOYDrwskTXfnHgUuUdV/zWtsq0TBMIxiOOean6VqW+jU6yhN20JWD6EJTJ8kLJm5T7rMT7NFVZfmuZ6qfgP4hoi8Aq994RSXdJCqPiIihwDfF5FfqOoDea5VqyiIyGzgC8DReC7Rfwd+CVwLDAHDwB+p6uOdzmVeQn+zbSRdbXBwoD8HNxZVu+/mPFHCEFfD90UmPP4hi3CEeyI2xFsoBFX9oYgcKiJzVXWLqj7i9j8oIjcDxwG5RKHuYPyngG+p6hHAC4B7geXA91T1MLwGluU12mf0AGkFIWvefuTqM48fW+LSg+QdvBYs7MucNiM8q3EvISJLRETc+vHAVGCriOwvItPc/rnAS4F78l6vNlEQkf2AVwBXAqjqLlV9AjgdWOmyrQTe2OlcUydLWWYahlEAwUbnju0AOYSmRwv9rwI/Bn5PRDaIyLki8m4RebfL8mbgbhG5E/gM8N9cw/PvA2tEZC3wA7w2hdyiUGf46BDgMeCLIvIC4A7g/cACVd0EoKqbRGR+1MEich5wHsCCAxaP7Z+/b8qZC41W0E3Nf9vIM30bRgoS7taZdgDY1WcenzmM1CkcVOfkenWjqm/rkH4pcGnE/luB5xdtT52iMAU4HvgzVf2piHyKDKEi13q/AuCIY45VMEEwjLR0KwhRpJ3mIk5IihKEOC/BBrNmo842hQ3ABlX9qdu+Hk8kNovIQgD399Ga7DMaTp72gX5rWzhi4cC4JS9FTXMRFoSiJ9gzQchObZ6Cqv5WRB4Wkd9T1V8CJ+M1ktwDnA1c4v7ekOZ85iU0m/WPpasNLpmXXCjUVZiPblg3tr5nYEFsvkkjm8dtT1l8VGk2FUUR8wYVMSleWkGIErUoL8EEoTvqHqfwZ8BXRGQq8CDwTjzv5ToRORd4CHhrp5PMnDqlY2Fi9D5NFwQ/PSwMTaCMCePCbQtJwjC8eWRCoZ9FRJK8mzSCYN3U01OrKKjqnUDUoI6Tq7bFMDrRSRCC+XxhGN2wrie8haJI6zFU1bA8f9+pPPr0rkqu1Rbq9hSMPmFw+mS27djd9fF1eQm7tm+BlGIQpKkeQ9FE9UTqJAxZBaFbL8EPKZfhJUyZOpk5z5uVLnOPPQZ1D14z+ojB6ZNLPf+2HbvHLU0hGH5qK+GCO6p9YHjzSGGCEDcQLSwIFlbOjomCUSndCEOSl5AkAt0IxK7tW8YteUgbbup1/J5IUcKQtvE4Kl+SIBjlYeEjo1Gsf2xkXO2ukyB0IkmEqhzE1va2heWvOqzQkcomCPVhnoJRCcHCt4gwUl5BiCKvZ9DvxE1al8VjgHyC8ND2HeO203aFNvZiomDUQlKB7b/I3TYuD06fnFoQym7ADoaQ2t62UOZ3m7N6CNbjqHtMFIzKCIdq4grvTgV6nJeQRQzCZPESnppkIYwyiBttnVUQfG/BF4Zbh7fmN66PMFEwaieqIM9ag88TkuokCE9NmjFuyYp5C50pqw3BPIbsWEOzUSmDAzMiC3x/HENWL6GI9omnJs1g5p5nJuwri6AwtLnxOYk08y/lEYSHtu8YN6r51uGtNqo5JSYKRmNIU8CXPdbBxxeJoDiEhcPojk6C0As9jCZPm8z+h8xOl/m2cm0pGhMFo+9JKuyLEIJ+GNnsc+Llt7Ds2AOA7N9eTisGURPdhXsdGd1jomAYJdFPYuDTqetpNw3JaWY7tRlRi8NEwehrygoJ9aMgnHj5LeNE4YiFA+O8haAgpPEKrKCvB+t9ZBgF04+CcOxf3Ti2HpyXyBeCtB/2OWjW9LHFqAfzFIy+oszG4n4UgyBDC6LHGUyYE6mGz2Zaz6P0mCgYraWq3kL9LgYnXn4Ls+ftmyqvfSGt+ZgoGLUSVXBnHSNQV1fRfhcDH78dwfcIhrel+39UJQZN9xJEZBnwKWAy8AVVvSSU/jzgKmAesA14u6pucGlnAx9zWS9S1ZV57TFRMBpHk8cDmBCMJ/yBnSSCXoJ5Bx4iMhn4DPBqYANwu4isVtV7Atn+HviSqq4UkZOATwLvEJFB4K/xvl6pwB3u2Mfz2GSiYNRGkwt/HxOBZNI2IKcdg+B/HCeOrNNWNN1LAE4A1qvqgwAisgo4HQiKwpHAX7j1HwD/6tZfA3xHVbe5Y78DLAO+mseg2kXBKeUaYKOqvl5EDgZWAYPAz4B3qKpNYNIymiYIVvino9M3E6IIC0LYS+gkBOG8TZjPaPK0fZg1lPojSnNFZE1ge4WqrnDri4CHA2kbgD8IHb8WeDNeiOlNwICIzIk5dlFao+KoXRSA9wP3Avu57UuBy1V1lYhcAZwLfK4u44xi2TbyzDhB6FQY5/16mRX25dDNVBR5xCB8XBOEIQNbVHVpTJpE7NPQ9oeAT4vIOcAPgY3AaMpjM1OrKIjIYuB1wMXAB0VEgJOAM12WlcAFmCi0Bl8Q0hbWPVWob90YnzYndwWudi75wf21ikH4HD0mDHFsAA4MbC8GHglmUNVHgP8KICIzgTer6nYR2QCcGDr25rwG1e0p/CPwYcAPTM4BnlDVUbcd6w6JyHnAeQAHHXRQyWYaRbBt5Blm0mMFfRqSxCCYp8eFoVMYqBNFiEH4fC0QhtuBw1zYfCNwBnsrxQCIyFxgm6ruAc7H64kE8G3gb0Rkf7d9qkvPRW2iICKvBx5V1TtE5ER/d0TWSHfIxeRWACxdujS3y2SUz8w9z4wXhDSFKTSvME1rd4tYtXbvby5CDOJmu03zmdXwuXtZGFR1VETeh1fATwauUtV1InIhsEZVV+N5A58UEcULH73XHbtNRD6BJywAF/qNznmo01N4KfAGETkNmI7XpvCPwGwRmeK8hQmulNGbjPuQTdZCNWv+okWkKBHocW8hbxgo69ToWQWiV1HVG4EbQ/s+Hli/Hrg+5tir2Os5FEJtoqCq5+NcHecpfEhVzxKRrwFvweuBdDZwQ102GsWSFDYa3TQ8tj5l4VC+CzW5Jh8ShtEN6xr5oZ2gZwDjBaEMMTCaQ91tClF8BFglIhcBPweurNkeIydxXkJQCIIUKhBGbnxBCIuBFfbtpBGioKo341rN3SCOE+q0xyiWsIcQJwZRxOXtabHoEW8BJgqCCUH7aYQoGO1l3Ifqt26MLeR3bnwIgGmL0vUkSxKWnhaMBhInCFV+19qoDhMFozTCghDEF4Ewwf1pBWLCdTcNT9jXdKHYtX0LU2fNrduMcRw0azrz9506TgyK/HRpL4vI5Gn7MGtJ73YaSMJEwSiUMSGIaOz1C+s4QQgTztetSASvHaRJQjFpZDM0TBTCFD01ycw9zyQKw+D0yR17IPVyd9SmYqJgFEaRghBF1hBTJ2pt0I7ontokbyHsJYQFIdhOlHcqkjj6pUtq0zBRMIqjQ8+ioCBsXz9RONK640WEmMKkbfwuUzya4i3cOrx1XE+jTnNVTRrZ3JUw9HL4qM2YKBiFMLr2pr3rm4bHpXUSg3Ballht0d5DJwofTxHyFprUE2lw+uTUkxd2KwxG8zBRMAojWGCGw0RJYhBm+/qNmRvxqhYH2Pt7i/Ye6hQG30sI9zRKM19V1cJg7QnlMKluA4zeZ3TtTbGCsH39xgmCsH14M9uHkwuZLCISZOfGh3K1W3RDlnEXE2jY6OuodoQsExi2brLDPsREweiK0bU37V02DY/t7xQqCopBWcIQtqMKcglD1PmC3XkrYv1jI5Vdq1N7Qt5G5vn7Tu2Fr641EgsfGZnx2w+6bTvwSfPlqm5CSUF7qgwndU0DJ8pr2pfx0mJCkB/zFIyuCHsHviBEhYuCZPiEYSFU7TEUSR3eQhXTWBTV66jo7zMYHuYpGJnY8a0V47azegc+VYpDz3gMLaDOHkhVegmTpk5t7TNlomB0TZIgxLUXzBpakFkQ8oSQfKronTS6abj7nkgxIaQmdVEtgjReQjftCRY2Kg4ThYoJhgR67WUPeglxgtCx8Xh4c+UhpCBt9xq2jaRrCxgciC+ce6k9ode/vNZErE2hQuqIERdFGkEwyqGM5yZKPMpuT7ARzL2BiUKN9JpIpBkD0MkL6NZLMPEplyQPI207gY1obgcmCjXSKy/R9is/FikGcQV1XMFfZ9goSC/2SKrKW4ij07Na1LNsk+DVj4lCzYz7VGWLaYogGOPZNvLMOHHopj0hjSCUGTrq9UZmEVkmIr8UkfUisjwh31tEREVkqdseEpHficidbrmiCHusobkmgi9Sk6ZMDrP9yo9F7+8QzilDBIrohdTL1N0Tac/AgkqmzM7K+sdGWDJvoG4zukJEJgOfAV4NbABuF5HVqnpPKN8A8OfAT0OneEBVjy3SptpEQUQOBL4EPBfYA6xQ1U+JyCBwLTAEDAN/pKqP12VnPxInBEb9hL9ZMeUFp1Z6fZsiu3BOANa7b9MjIquA04F7Qvk+Afwt8KGyDaozfDQK/E9V/X3gxcB7ReRIYDnwPVU9DPie224d4YnDmhJG6hdBaGS7QtrJ8RLyDQ7MiF060UtdUXuMuSKyJrCcF0hbBDwc2N7g9o0hIscBB6rqv0ec+2AR+bmI3CIiLy/C2No8BVXdBGxy6yMici/ezTgdONFlWwncDHykBhP7kiaHaJpqV6FknAdpdO1NhXgLJggZ2WdaloGKW1R1aUyaROzTsUSRScDlwDkR+TYBB6nqVhF5IfCvInKUqj6Z1rAoGtHQLCJDwHF48bIFTjB84Zgfc8x5vvI+9thjVZlaKE3zFh46/521Xt9wJHkMOabaTuMt1E0Vcy81jA3AgYHtxcAjge0B4GjgZhEZxouqrBaRpaq6U1W3AqjqHcADwOF5DapdFERkJvAvwAeyKJyqrlDVpaq6dN68eeUZWDI2/3w6yvASGhlC8tm6MXqJIPjVu24wL6FWbgcOE5GDRWQqcAaw2k9U1e2qOldVh1R1CPgJ8AZVXSMi81xDNSJyCHAY8GBeg2rtfSQi++AJwldU9etu92YRWaiqm0RkIfBofRZWT9aeSHEFQtaQQpO9hL4IG1XE4MCM3F1QjeJQ1VEReR/wbWAycJWqrhORC4E1qro64fBXABeKyCiwG3i3qm7La1OdvY8EuBK4V1UvCyStBs4GLnF/b6jBvMJJGnw0wVtIKQpJNcSiYs1Geor+LGcW2vT/Hpw+ua8GsanqjcCNoX0fj8l7YmD9X/Aq1YVSZ/jopcA7gJMCgy9OwxODV4vI/Xh9dy+p0cZaaOr0F7POvWjvesm191lLFo0tRrH0QtuCUR919j76EdEt7wAnV2lL2ZRRyOeNIxfBrCWLCp2TqGoBKHK21Dq9BJ+ivYXg+II8YaaZe57JPFah37yFJlF7Q7MRTdO8Bd9LCHoLuc4X8ATq8AjaPH12GgYHZsQW9E9NmjGhEM87AM3aLnoHm+aiZNrqJfjEeQtNDvu0WRCyTIUxddbcCd2gkwr/pybFC0kasnoM5i3Ug4lCCRQlBHXPdeMT9g5mnXvRuJHPTRaAMGUIQhNCR2Ns3choTFKnZylNgW3C0H5MFIy+oM3eQVqiKhnB7s+DCccGu7H6hXpVISF/QFuUOAxOn1xPw/nkfTKNPO8lrE3BmECaxuOi2hbKZtqig0wQAhQZzvTbHoJLGroVk8HpkycskO27EEZnzFMwWokJQbGEB73FkbbHUjc9koxqME/B6JomegtVewaNak/w6TA/UlU92zp5D9YjqZmYp2C0AvMMyiett5AFXxjMa2gO5in0GFV1R007KG3WuRcVOoCtG0wQslPlOJg0Bb55Dc3BPIWGU2e31O3rN3LQJ7+YnOfKj9XaJbVuQRjdNNzMEFJJJPX0yetFWDtDMzBPoYdo0qC1JlC3IDSaFN9dqPJ5Kruwtx5IxWGi0APUNeVFGi+hLkwQcpLjYz1xFDFeIE8YyYShGEwUeoGtGxvnJZgg9DABQSj6uUoShrLHMRjF0Io2hdE9e+o2oXVEeQl1CoGPCUJ/0G37QmXewuR92DOwoJprVUwrRKFJFB7qyeHmZ5lKuVOoyDDSEvYWLKzTW7QmfGQPXn/QpO8qT1k41NyeRz0+L4/1QqqP1ohCE2iSl+DTtLaIttBoMegkCBULhu85pG0rMEGol1hREJElIvLSiP0vF5FDyzWrO2r3Fkro0WFMpG5voZGCkEYMGk6WSfWM8kjyFP4RGInY/zuXZgQYq5Fv3ViMOBQoMOYtFEdjBaGKY3LQqVeSiUFzSBKFIVW9K7xTVdcAQ6VZ5BCRZSLySxFZLyLL0x63beSZcUsVRBa6vjh0U7ibx9GROryFxglCXu+gAZ6FiUHnsk5EponItS79pyIyFEg73+3/pYi8pgh7kkRhekLac4q4eBwiMhn4DPBa4EjgbSJyZFz+3U3vkZqlkDdBSE3VwjC6abjS6yVSVIFeoTAEP+hTF037ilvKsu5c4HFVXQJcDlzqjj0SOAM4ClgGfNadLxdJonC7iPxJxI84F7gj74U7cAKwXlUfVNVdwCrg9KQD6vpnpw7NdCrsiwo7lUjcVNl1TqFdd/tCLTSght9rbNuxu3GC4EhT1p0OrHTr1wMni4hsk5KiAAAf7UlEQVS4/atUdaeq/hpY786XiyRR+ADwThG5WUT+wS23AO8C3p/3wh1YBDwc2N7g9o0hIueJyBoRWbNtq/fx8Qb/4z2CIaXw0iNEfa85an+VVCkMjfIWepCwt1DF6OUGlAlz/bLKLecF0jqWdcE8qjoKbAfmpDw2M0mD174InIXXfnC02/cfqvr9vBdNgUTs03EbqiuAFQDPP/b4cWn+Q+B/rq8srAG3Oezc+JCNdu4Rps6ay67tWyq5VlmCsFsmZWkP2aKqS2PSOpZ1CXnSHJuZJE/hauDbwEuAK1T1nyoSBPAU78DA9mLgkawnaUANoZXEeQd1f4mtKo+hVm+hhaGjsryFHnn/05R1Y3lEZAowC9iW8tjMxHoKqnqdiPwH8HFgjYj8M7AnkH5Z3osncDtwmIgcDGzEa0w5s8TrGRnY8a0VgPuWQk5hKHo+pao8Bl8YKu+RtHVjK4RhQhgp5Dn0Ua+kNGXdauBs4MfAW4Dvq6qKyGrgGhG5DDgAOAy4La9BnUY0Pws8DUwDBkJLabi42fvwPJV7getUtZ75o1tCU0NdZXgXVbcxVO459FAblJFMXFknIheKyBtctiuBOSKyHvggsNwduw64DrgH+BbwXlXN7R7Fegoisgy4DE+ljlfVSocLq+qNwI1VXrPtZJkgLw7fS/Br4zu+tYLpy85LOqQWqm5jqNxzaInH4BNuZ+inr7BFlXWq+vHA+g7grTHHXgxcXKQ9SZ7CR4G3quryqgXB6C/KaovYufGhWsYyVOY5mMeQSNkdTdpKrCio6sstZNM+8oSR4o71vYem0npxaAl1dFc1JmKzpPYho2tvGluaQhU9l+oY6Fa6OBQxk25Nn3s1momJQp+TVhjKFpCqvupW1who8xyMXsG+vGYYFTK6abh5E+vheQtTFh9VtxkTaGqD8+49PTMOIjPmKRgdvQA/vayabtXffq57viTzGLJhbQvVYqLQRPwpkVvU5dAYTxOFwdoWDDBRKIRC48VhIahIGOK8hbZ5CT51ewvQTGFoKuYtVIe1KRRIKfHiCvuiV9UbqS4haCKFPDMtG8hm1IuJQgFMWThUXK2vgQOSon5bGaOY/RHIVdXibWbVZlHVzKlGMhY+Kgi/kLSQQHcEC+d+K6jtmTGahIlCAeSdT6jJVOElRIlAVcLQhLaFQmigh5kF8xKag4lCgbTNWyjrd0xbdNC4xTDKwuY/yo6JQk7CXkITZwzthjhBqPL39ZO3UIgA5/QW9gwsqKXGnvaa1gOpGqyhuUuS+nQ3ccRqVgptPDeqw3oiTWBw+uTCRx8/u2cPjz69q9BzNgXzFPLQZy9fUV5CW7ypoqm7B9ukkc1As+P75i2Uj4lCF3Q78tMf5NYrk6M1wePppxASFCwMPd74HEc3wmBtC+lptSiU4d5lrUUliUCviAMUX7s3byGeQp+JLoWhyd5Ct5gwpKPVolA3aV/uJovDlIVDjfAYjBy00GOwMFJ51CIKIvJ3InKfiNwlIt8QkdmBtPNFZL2I/FJEXlOHfXGUXXtqrDCUNA4jrbcQF0IquktrU0JI0IxnoenegoWRyqEuT+E7wNGqegzwK+B8ABE5EjgDOApYBnxWRHL9F28d3prT1GppQmEQpOyBeVmEIW7Uc1Hi0OoxEy30FsATBvMaiqUWUVDVm1R11G3+BFjs1k8HVqnqTlX9NbAeOKEOG8P4tSa/h0aZ9FJjdNUkCUA34tDkQXR1/P/Dz3fZ3sK2kfEFerfvV1vFQUQGReQ7InK/+7t/Qt79RGSjiHw6sO9mF3W50y3zO12zCW0K/x34pltfBDwcSNvg9k1ARM4TkTUismbb1uQHt9e8hSB1CkNV03eUMW1G2qWvaJi3EBYEn0kjm3OJQ8tYDnxPVQ8Dvue24/gEcEvE/rNU9Vi3PNrpgqWJgoh8V0TujlhOD+T5KDAKfMXfFXEqjTq/qq5Q1aWqunRwztzif0CAbryEIgtz8xiMOqmrbcG8BsCLnqx06yuBN0ZlEpEXAguA3PPflzaiWVVPSUoXkbOB1wMnq6pf8G8ADgxkWww80q0N8/edOtYt9dbhrbxkaE63p6okbNQkqp7kb/qy89jxrRWVXrMXaOo3nfOy/rGRCY2+Ue/YpJHN7BlY0NU1yvy+867dykPbd6TNPldE1gS2V6hq2od9gapuAlDVTVHhHxGZBPwD8A7g5IhzfFFEdgP/AlwUKG8jqWWaCxFZBnwEeKWqBiV9NXCNiFwGHAAcBtzW6XyTmxAEc1T1oZpcRI3EblhowWg3viCkqdHnEYaGsEVVl8Ylish3gedGJH005fnfA9yoqg+LTAi2nKWqG0VkAE8U3gF8Kelkdc199GlgGvAd9yN+oqrvVtV1InIdcA9eWOm9qlrspCVdsmdgQWRNZtf2LUydVV74qpBaYprpOOYsMmFoKIV6Cw2YGymuLSGJboWhTG+hKJKiKiKyWUQWOi9hIRDVJvCHwMtF5D3ATGCqiDylqstVdaO7xoiIXIPXcad5oqCqSxLSLgYuznrOMia9ykqRHz7PXQh08+LXLAw2CV9/EfQS0oRn2ywMCawGzgYucX9vCGdQ1bP8dRE5B1iqqstFZAowW1W3iMg+eOH673a6YIMCL80n7oEsuiGua0GYs2jv0i15jzdKoS1i2Y2X0OdcArxaRO4HXu22EZGlIvKFDsdOA74tIncBdwIbgf/b6YKtmjq7Tm9h1/YtYwqb5wXuShCsEDcKoqrYfVYvIZi3n7wFVd1KROOxqq4B3hWx/2rgarf+NPDCrNdslSjUwbgHOmfoxQTBaDO+l5C3u2i/CUPVWPgoA2lqNJW6+S0UhDZ2vyyCtt6Xbrp693hPpMbTOk+hthDS1o09JQjhF6vfxmH0CoWLQU0ViaK8hLzUff1eoHWiECY4gK0w6u66mePFjqtlxXW5NaqnLq+g6hq4PW/NpPWiUBRjD3CMIIS9hOA0zGnm2Cm7IEjzwpsw1EtbQ0RRo5e7pSmho12jexje1k6vw0QhCyFBiAoXlT4nfxdeQpYXqSkvXRqKGNeQpyAuKlzYVjEAb3qZ+ftOnbDfKh/NpZWiUGS7wraRZ5hZyJnqoZcK+Sz4BWmwQM1aSOctjPNcu4jrZ6aG9oQoQeiWtj7LTaOVohDGb1fIOylekLYMJoqj7Ok7uiWpIPXT0vxvii6QG1/br7GnWlToKG2o0oSgeqxLagL9Pvqy8imTXcEVVcBm+VZ0p3yNL8CLJMMI9bhCuqjKQbjnT6cC3wShHlorCuHaie/GFvHBnbZ7CU0iixgUcVxr6PHpSkwQ6qMvwkc+WeKbRfaYKJSKZ7msK4yUqUD370dEz7C+EoYSnou8//tO71A4jGRiUD99JQo+6x8bARgbv1BkY9i0RQeV3wOpYnZt3zJueoBwGOCpSTPG7atMRMKFYII49CwV1/br+HaBCUGz6EtRKJus3/7tqlZcMVnmjenGuxibdjzt9N1J96GXvw1Rc8gnWEAXLe7hykPV2LxH6TBR6IIsvVw6nSMTGQuMMmpgUS91pS96zD2Ysvio8d+zqFMYejiWXyRRHTWKEAYr3MultaKQZpxC5ikwQgVNN+JQ1UyoVQlCHFm8heC044l0uA+RwgDFiUPLC/uinxn/HQy3K2QRhqYKwI5nd3PfppG6zSiF1opCWtK2JyT1qy60MbOAgqeVMdqE+zJl8VFj65H/pzxeQ8uFoCz8djvwxCFKGIxm0veikJkia58lFDhNEoQ03kKqsRAZ7tPUWXPxfb9x4pD1/9ZnYlD2cxMlDEYzqXWcgoh8SERUROa6bRGR/yMi60XkLhE5vk77EmnxJy+fmjSj3ppc8N6kCBnFEVnQdbr3Df/fZCFtQV+0IAS9hCDbduyu/TvqRmdq8xRE5EC8b44G+2++FjjMLX8AfM79TWT3nnRtCFE8+vSu2BCSX7vxY6DBl2fSyObCCo8ivm2Q58WOE4AswhAXI841zqHAUFpkSKnH6OZ/XLXnGCcIQYLvajfegwlLudQZProc+DBwQ2Df6cCXVFWBn4jIbBFZqKqbyjDAb2ROEgYfv4D0C78yX7ZO5y6iL3nRnkD4/gTpShgKLrSbOC14k0J9dWEFfPOoRRRE5A3ARlVdKyLBpEXAw4HtDW7fBFEQkfOA8wAOWHxgZhvS9joK96CIKky76WKXVCh3Ol8ZXkFRJA1yq5tYr6HEa/UTcV5CGYNE+wURGQSuBYaAYeCPVPXxiHyXAq9zm59Q1Wvd/oOBVcAg8DPgHaqaWPiVJgoi8l3guRFJHwX+Cjg16rCIfRp1flVdAawAeP6xx0fmyUKnBzdcowm6vWXVun3a0K+7SbOu5vUa+qHAnzSyGQr4fwUrX+GKmIlEKpYD31PVS0Rkudv+SDCDiLwOOB44FpgG3CIi31TVJ4FLgctVdZWIXAGcixeWj6U0UVDVU6L2i8jzgYMB30tYDPxMRE7A8wyC1f7FwCNF2tXJQ0j74Jbp9sZ138siDk0QgrpHsCYRVbBHCUU/CECYtIIZHpwW/o5J1nctSJJgFP553WZzOnCiW18J3ExIFIAjgVtUdRQYFZG1wDIR+RpwEnBm4PgLqEsU4lDVXwDz/W0RGQaWquoWEVkNvE9EVuE1MG9P057w7J49pT0oWc5bVM0nriGuCQV9HsLeQhlTc0edM63n1RYBiHtO0gp0t55UFkHoRMsK/rkisiawvcJFOtKwwC8DVXWTiMyPyLMW+GsRuQyYAbwKuAeYAzzhxAL2huMTado4hRuB04D1wDPAO+s1JxtFPsi+wOTtqVE3ab2FbgqipO6oRdhUFVWJfZb2nk4CGfYSkgThoe07Es910KzpielNZOezuxnenHpE8xZVXRqX2CHU3hFVvUlEXgTcCjwG/BgYJUM4PkjtoqCqQ4F1Bd5bnzWd6fSAQ7qHPO48/rHBF6tbgegmxFWG8Iwr9EaeYXBgRu7PnEYVWmm8hPD+ssWhqd5dke1WcYKQ5l1Jky/ufUp7/qYTF2oHEJHNfg9MEVkIPBpzjouBi90x1wD3A1uA2SIyxXkLqcLxtYtCEezarYkPSJ6aSDcPXp6HNXhsWCCC4amy2jTizluWWAQLo6TeQWERKKrRuqmFdtVE3YfBmLxBLyHueSmywG5L4d8lq4GzgUvc3xvCGURkMjBbVbeKyDHAMcBNqqoi8gPgLXg9kCKPD9MKUehEGQ/V8LZya5hDgzPG7E4Sh6pIEqGsghEsVKJCOX4oKUuMP4uXYJSD/3yG37es78rQoP3fAlwCXCci5+IN9H0rgIgsBd6tqu8C9gH+03XceRJ4e6Ad4SPAKhG5CPg5cGWnC7ZCFHaN7qnsWlEPeHC2xCMWDhR6nSRxgHiBiGvfKENQkrrrdkuSIDSla2s/EuclRAlCtxWnpOOCglF2xawJqOpW4OSI/WuAd7n1HXg9kKKOfxA4Ics1WyEKEP+AFFXr6CQGSfu6wReX4W3PjP2GsDhA9sbtbrsBZiFr20fWht80vZiCXkLWUFtRobKmjdbN8ru8/0fyuxPsgupPQ3/QrOmJnnnc+5G2MtUPQlA3rRGFONLWOrIcV8U86vdtGokUBohudyiCMgQjrUDkFYbwufLQtMK8KNL+L8pqfE96b6LSivK6jWy0QhTCH7wos9YR9fB26po2tKC7hzssDDBRyPK2l6QVlSIEo1OYKYsw+CNuRzesYxJeuClKDNpawOclbipr//6nHdEc5S1E0U1Fqq0fsWk6rRCFMEXXOuIezrT9lMP5kkQimHdowcA4YYCJQpY3PBYlKlm9j27bL6K+zJWlq6j/lbUsDdJZw21Ft8EUMZalyDBf8N6PE4StGyHlOJDwSOY0xL073VagjOJohSjsfLbzA5k1lplUS8kwaKXj8cGXIHze4c0jY8LgE7a32xhrkpgUIRSQrkEc4r/MVfS0Ht0UyE0cWVtkmC9874uYLDDcrhB8fju9O1kqUJ3I+572K60QBei+5pHVRY26zhOPPd3xuNnz9k19vnB68DcU5QVl9TjyCkUngUj6ZGOSOCR5CUVOu9ArRP3OLEIRHFk+uvYmprxg/LyVeb6g1k0h3dSCfffonlTvfS/SGlGII65Wnuc8QdI+GE889nSsMKS5dpLtWYUtSkSiPI5uhCJInGjETTqYds6n8AeP/PSoEEZUIVlVO0wn0tqR93qdJnkMf0wKGPuexeiGdWPTifhTY4eFIXzf0/6uuHen2/fEKIZWiEJQtZMeqCzeRKcaSje1hLzCkEQWwUvrbeTt5ttpKg+fTiO2kzyItILQyyNsiwrn+cTd7zFhSDEPVd4G/KT3J5hWhEC0tUZfFq0QhSDhByDNQ5XVRc3zkOURhiTS/oY48QgLRVJIqttuvj7hQq7ToLxOHkS4gEo7IVvRjfZVUsS0LmFxCAuDz7aRZ7oWAf+5Cj6fWd4fK9Crp3WiEKYbkUh7rjznSWtHp2vOnrdvKi/JJ0o8ooSi27aLNA3fcd1q4wblJQlEltk5O9lW1sCoTmLTTeguibQemk/wc7RBYSiDqOf5iUfHP2uz51sPpDpphSjsHt0z7sFKeqiSCtlwoZpGBMIP9IRzxtiSRhhSXT+QJ42ARJG2x0fe0ahj1wsUgnGD8qDznE9pZ+esexRsN9cvWiggWnx9wsLg03WjcsD+JC826v1J+y5npdO7ani0QhTCdPtQZXJrUz5gTzw6kksYiiSt15S1W2CWhu5OXWqTRm6n6UGUZu6dvIOi6hppW9RULlnCd0lho7hGfN/O8H0OP39p3qGoPGneaROA7mmlKAQpo9aR9YHrRhg6CVQaGzr93m5FAvKN0g6TZnBe0tQeaTyDIkfH5jlXlKDknVCxm5BdkKTp2jsRd+/v2zQS246Qt8C2Ar9cWi8KQfLGLvM8jJ2EAfYWykUIQlS+okQCuus/niUsFZ7eA6IFYpxNocKxidMkdLKpU3q3XkpaL6ObnlXBc4cbltOIwcgj68dtDxywJLMNRnG0QhR279ox7sFK+1BVXeNIEgYopg0jy7FFikQasnSrjeoNlTY2n3dakqIoY8qGJNEo0ssIi0Wnex81atl/ftKKQXh/0eIQdz1jPK0QhTDdCERVdBKGTscWbUuQrCIRpgzR8AvWpC6zZU5Jkoes184rImk9ozy9yNLc6yTvIEvBHJU36/tclhCEO7e0idpEQUT+DHgf3gem/0NVP+z2nw+cC+wG/lxVv53nOk10TbsRhqy1rTBpfne3jXpjxxfUZTcoLnGN3kUIQZS93fYIK6LDQBq7i/A+igirdZruJfgsFVkwW22/fGoRBRF5FXA6cIyq7hSR+W7/kcAZwFHAAcB3ReRwVS1s/uMyXdMs58zjMXRDt7+7rO6BiddMGNHabc0/y3QkRZ6/6N5lZU3TnueaSb2KrBDvPeryFP4UuERVdwKo6qNu/+nAKrf/1yKyHu9Tcj8u2oC84hD1sGcNW6UVhiLd1DyhtTg7yhSLbqc8aMpI2G68kTxUES6Lu7d5QkVGc6hLFA4HXi4iFwM7gA+p6u3AIuAngXwb3L4JiMh5wHkAk6bP6tqQNHHLbh7utKLjv0ixPZNKjFsW1fZStI2demn1Oll/R50TxHXbE84EoXcpTRRE5LvAcyOSPuquuz/wYuBFwHUicgggEfk16vyqugJYATBl1qLIPN1SRgw0rTjURZPaXvK2baQ5X5kU7TmlFZFuxCOrQHW6lyYGvU9poqCqp8SlicifAl9XVQVuE5E9wFw8z+DAQNbFwCOdrrXn2Z2MbHoAgIGFh+YxuzSytjnUTRE9P4okzXQidQurT9Gilvq6JXhSSfc0jwBU+b761+pFROStwAXA7wMnqOqaiDzTgR8C0/DK9OtV9a9d2tXAK4HtLvs5qnpn0jXrCh/9K3AScLOIHA5MBbYAq4FrROQyvIbmw4Dbspx4ZNMDJgwlUYQ3UVZ34aYIQhx1tMd0SxZvYOuNHx1bn3PaxdH5Q+9ksJAuSxx6WQhC3A38V+DzCXl2Aiep6lMisg/wIxH5pqr6ofi/VNXr016wLlG4CrhKRO4GdgFnO69hnYhcB9yD11X1vUX2PGoCRQvDwAFLanPZ04TGkmxrUsiq0z0sy7ZuJ1Ss0oYg3Txr4QJ/YOGhpRbaLRIEVPVeAJGoyPpYHgWecpv7uKXrkHotoqCqu4C3x6RdDERXOYxGUpQolRGyarJtaWiSBxS+B0Evwd+O8xbC+EJRdAHeUEGYKyLBsM8K1yZaGCIyGbgDWAJ8RlV/Gki+WEQ+DnwPWO73+oyjlSOam04Z3oJ/3rbR5N/U5JHzRdNJELqlqLBR1WIQnlqnA1tUdWlcYlKnHFW9Ic0FXETlWBGZDXxDRI5W1buB84Hf4oXoVwAfAS5MOlcrRaHOdoW0BXQZg+jaLA5Np0mhsKLo5jnK4i1EXjOicI97lxvqFWQmqVNOF+d6QkRuBpYBd6vqJpe0U0S+CHyo0zlaIQqT9pnWuMbltLH+umLZRrn0mkhkEYCivIQosrzH3b7z2xL73vQeIjIPeNYJwnOAU4BLXdpCVd0kXqPEG/EarhNphSjURZkvRz+Rp2bZKzQ11FS0V9kP/8sqEZE3Af8EzAP+Q0TuVNXXiMgBwBdU9TRgIbDStStMAq5T1X93p/iKEw0B7gTe3emaJgqGUTFNGAPSrRhYRahaVPUbwDci9j8CnObW7wKOizn+pKzXNFHoEns5iiNvHLoNdNXVsyHTSPv0+/+wLZgoGEaPUnWHAqsI9Qd9IQppa1RpXzJ7OYrHvIXexv537aEVonDsYQtZYwV1z5NGbK3wiSbNdBNZz9PrxP0WkY9VbElv0QpRMPoH8yiMJhCchLNtTKrbAMMw8hE13cTAAUvGlm7P08u06bdUjYmC0XPYC7+Xou6F3VPDx0TBMFrI8Bf+29h6kwbLVYEJXD5MFIyepMkvftawTbeYl2CUgTU0Gz2LFWbJBL2FsrH/RXswT8EwDMMYw0TBMAzDGMNEwTAMwxjDRMEwDMMYQ7xvPvc2IvIY8Ju67QgwF9hStxERmF3paaJN0Ey7mmgTxNv1PFWdl+fEIvItd/40bFHVZXmuVyWtEIWmISJrkr7JWhdmV3qaaBM0064m2gTNtavpWPjIMAzDGMNEwTAMwxjDRKEcVtRtQAxmV3qaaBM0064m2gTNtavRWJuCYRiGMYZ5CoZhGMYYJgqGYRjGGCYKBSAiwyLyCxG5U0TWuH2DIvIdEbnf/d2/AjuuEpFHReTuwL5IO8Tj/4jIehG5S0SOr9CmC0Rko7tfd4rIaYG0851NvxSR15Rk04Ei8gMRuVdE1onI+93+uu9VnF1136/pInKbiKx1dv1vt/9gEfmpu1/XishUt3+a217v0ocqtOlqEfl14F4d6/ZX8j9sBapqS84FGAbmhvb9LbDcrS8HLq3AjlcAxwN3d7IDOA34JiDAi4GfVmjTBcCHIvIeCawFpgEHAw8Ak0uwaSFwvFsfAH7lrl33vYqzq+77JcBMt74P8FN3H64DznD7rwD+1K2/B7jCrZ8BXFuhTVcDb4nIX8n/sA2LeQrlcTqw0q2vBN5Y9gVV9YfAtpR2nA58ST1+AswWkYUV2RTH6cAqVd2pqr8G1gMnlGDTJlX9mVsfAe4FFlH/vYqzK46q7peq6lNucx+3KHAScL3bH75f/n28HjhZRKQim+Ko5H/YBkwUikGBm0TkDhE5z+1boKqbwHvZgfk12RZnxyLg4UC+DSQXQEXzPufGXxUIrVVukwttHIdX02zMvQrZBTXfLxGZLCJ3Ao8C38HzSp5Q1dGIa4/Z5dK3A3PKtklV/Xt1sbtXl4vItLBNEfYaAUwUiuGlqno88FrgvSLyiroNSkFUza2q/smfAw4FjgU2Af9Qh00iMhP4F+ADqvpkUtaIfVXaVfv9UtXdqnossBjPG/n9hGtXYlfYJhE5GjgfOAJ4ETAIfKRKm9qAiUIBqOoj7u+jwDfwXprNvnvq/j5ak3lxdmwADgzkWww8UoVBqrrZvdB7gP/L3pBHZTaJyD54Be9XVPXrbnft9yrKribcLx9VfQK4GS8uP1tE/K83Bq89ZpdLn0X6EGIem5a5EJyq6k7gi9R4r3oVE4WciMi+IjLgrwOnAncDq4GzXbazgRvqsTDWjtXAH7teGS8Gtvuhk7IJxXLfhHe/fJvOcL1XDgYOA24r4foCXAncq6qXBZJqvVdxdjXgfs0Tkdlu/TnAKXjtHT8A3uKyhe+Xfx/fAnxfVQutlcfYdF9A1AWvjSN4r2p53nuOulu6e30BDsHrAbIWWAd81O2fA3wPuN/9HazAlq/ihReexasZnRtnB547/Rm82PAvgKUV2vTP7pp34b2sCwP5P+ps+iXw2pJsehle6OAu4E63nNaAexVnV9336xjg5+76dwMfDzz7t+E1cH8NmOb2T3fb6136IRXa9H13r+4GvszeHkqV/A/bsNg0F4ZhGMYYFj4yDMMwxjBRMAzDMMYwUTAMwzDGMFEwDMMwxjBRMAzDMMYwUTB6Hje76K9FZNBt7++2nycih4vIjW52zHtF5DoRWVC3zYbRVKxLqtEKROTDwBJVPU9EPo83c+3leH3SP6iq/+byvQp4TFXvjj2ZYfQxJgpGK3DTQ9wBXAX8Cd5kcm8HTlTVP67TNsPoJaZ0zmIYzUdVnxWRvwS+BZyqqrvcBGl31GyaYfQU1qZgtInX4k2pcXTdhhhGr2KiYLQC99nFV+PN3vkXbmK0dcALazXMMHoMEwWj53EzYn4O7/sDDwF/B/w9cA3wEhF5XSDvMhF5fj2WGkbzMVEw2sCfAA+p6nfc9mfxPrRyAvB64M/cx+XvAc6hvm9bGEbjsd5HhmEYxhjmKRiGYRhjmCgYhmEYY5goGIZhGGOYKBiGYRhjmCgYhmEYY5goGIZhGGOYKBiGYRhj/H/Co7R7gdGHjwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "surf_mask_c = ds.hFacC[0] > 0\n", "ds.Eta[0].where(surf_mask_c).plot.contourf(levels=30)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating the grid object\n", "\n", "Next we create a `Grid` object from the dataset.\n", "We need to tell xgcm that the `X` and `Y` axes are periodic.\n", "(The other axes will be assumed to be non-periodic.)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mInit signature:\u001b[0m \u001b[0mxgcm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGrid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcheck_dims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mperiodic\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdefault_shifts\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mface_connections\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcoords\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mDocstring:\u001b[0m \n", "An object with multiple :class:`xgcm.Axis` objects representing different\n", "independent axes.\n", "\u001b[0;31mInit docstring:\u001b[0m\n", "Create a new Grid object from an input dataset.\n", "\n", "Parameters\n", "----------\n", "ds : xarray.Dataset\n", " Contains the relevant grid information. Coordinate attributes\n", " should conform to Comodo conventions [1]_.\n", "check_dims : bool, optional\n", " Whether to check the compatibility of input data dimensions before\n", " performing grid operations.\n", "periodic : {True, False, list}\n", " Whether the grid is periodic (i.e. \"wrap-around\"). If a list is\n", " specified (e.g. ``['X', 'Y']``), the axis names in the list will be\n", " be periodic and any other axes founds will be assumed non-periodic.\n", "default_shifts : dict\n", " A dictionary of dictionaries specifying default grid position\n", " shifts (e.g. ``{'X': {'center': 'left', 'left': 'center'}}``)\n", "face_connections : dict\n", " Grid topology\n", "coords : dict, optional\n", " Excplicit specification of axis coordinates, e.g\n", " ``{'X': {'center': 'XC', 'left: 'XG'}}``.\n", " Each key should be the name of an axis. The value should be\n", " a dictionary mapping positions (e.g. ``'left'``) to names of\n", " coordinates in ``ds``.\n", "\n", "REFERENCES\n", "----------\n", ".. [1] Comodo Conventions http://pycomodo.forge.imag.fr/norm.html\n", "\u001b[0;31mFile:\u001b[0m ~/miniconda3/envs/pangeo/lib/python3.6/site-packages/xgcm/grid.py\n", "\u001b[0;31mType:\u001b[0m type\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xgcm.Grid?" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "OrderedDict([('standard_name', 'longitude'),\n", " ('long_name', 'longitude'),\n", " ('units', 'degrees_east'),\n", " ('coordinate', 'YC XC'),\n", " ('axis', 'X')])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds.XC.attrs" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Z Axis (not periodic):\n", " * center Z (15) --> left\n", " * left Zl (15) --> center\n", " * outer Zp1 (16) --> center\n", " * right Zu (15) --> center\n", "X Axis (periodic):\n", " * center XC (90) --> left\n", " * left XG (90) --> center\n", "1RHO Axis (not periodic):\n", " * center layer_1RHO_center (30) --> outer\n", " * outer layer_1RHO_bounds (31) --> center\n", " * inner layer_1RHO_interface (29) --> center\n", "Y Axis (periodic):\n", " * center YC (40) --> left\n", " * left YG (40) --> center\n", "T Axis (not periodic):\n", " * center time (1)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "grid = xgcm.Grid(ds, periodic=['X', 'Y'])\n", "grid" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mSignature:\u001b[0m \u001b[0mgrid\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mda\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mDocstring:\u001b[0m\n", "Interpolate neighboring points to the intermediate grid point along\n", "this axis.\n", "\n", "Parameters\n", "----------\n", "axis : str\n", " Name of the axis on which ot act\n", "da : xarray.DataArray\n", " The data on which to operate\n", "to : {'center', 'left', 'right', 'inner', 'outer'}\n", " The direction in which to shift the array. If not specified,\n", " default will be used.\n", "boundary : {None, 'fill', 'extend'}\n", " A flag indicating how to handle boundaries:\n", "\n", " * None: Do not apply any boundary conditions. Raise an error if\n", " boundary conditions are required for the operation.\n", " * 'fill': Set values outside the array boundary to fill_value\n", " (i.e. a Neumann boundary condition.)\n", " * 'extend': Set values outside the array to the nearest array\n", " value. (i.e. a limited form of Dirichlet boundary condition.)\n", "\n", "Returns\n", "-------\n", "da_i : xarray.DataArray\n", " The interpolated data\n", "\u001b[0;31mFile:\u001b[0m ~/miniconda3/envs/pangeo/lib/python3.6/site-packages/xgcm/grid.py\n", "\u001b[0;31mType:\u001b[0m method\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "grid.interp?" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('time', 'Z', 'YC', 'XC')" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds.THETA.dims" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "array([[[[ 0. , ..., 0. ],\n", " ...,\n", " [-0.113758, ..., -0.320469]],\n", "\n", " ...,\n", "\n", " [[ 0. , ..., 0. ],\n", " ...,\n", " [ 0. , ..., 0. ]]]], dtype=float32)\n", "Coordinates:\n", " * time (time) timedelta64[ns] 11:00:00\n", " * Z (Z) float32 -25.0 -85.0 -170.0 -290.0 -455.0 -670.0 -935.0 ...\n", " * YC (YC) float32 -78.0 -74.0 -70.0 -66.0 -62.0 -58.0 -54.0 -50.0 ...\n", " * XG (XG) float32 0.0 4.0 8.0 12.0 16.0 20.0 24.0 28.0 32.0 36.0 ..." ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "theta_x_interp = grid.interp(ds.THETA, 'X')\n", "theta_x_interp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "theta_z_interp = grid.interp(ds.THETA, 'Z', boundary='extend')\n", "theta_z_interp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds.THETA[0,5].plot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "theta_x_interp[0,5,].plot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds.THETA.sel(XC=200, YC=45,\n", " method='nearest').plot.line(y='Z', marker='o',\n", " label='original')\n", "theta_z_interp.sel(XC=200, YC=45,\n", " method='nearest').plot.line(y='Zl', marker='o',\n", " label='interpolated')\n", "# plt.ylim([-500, 0])\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Kinetic Energy\n", "\n", "Finally, we plot the kinetic energy $1/2 (u^2 + v^2)$ by interpoloting both quantities the cell center point." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# an example of calculating kinetic energy\n", "ke = 0.5*(grid.interp((ds.U*ds.hFacW)**2, 'X') + \n", " grid.interp((ds.V*ds.hFacS)**2, 'Y'))\n", "print(ke)\n", "ke[0,0].plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Vorticity and Strain\n", "Here we compute more dervied quantities from the velocity field.\n", "\n", "The vertical component of the vorticity is a fundamental quantity of interest in ocean circulation theory. It is defined as\n", "\n", "$$ \\zeta = - \\frac{\\partial u}{\\partial y} + \\frac{\\partial v}{\\partial x} \\ . $$\n", "\n", "On the c-grid, a finite-volume representation is given by\n", "\n", "$$ \\zeta = (- \\delta_j \\Delta x_c u + \\delta_i \\Delta y_c v ) / A_\\zeta \\ . $$\n", "\n", "In xgcm, we calculate this quanity as" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "zeta = (-grid.diff(ds.U * ds.dxC, 'Y') +\n", " grid.diff(ds.V * ds.dyC, 'X'))/ds.rAz\n", "zeta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "...which we can see is located at the `YG, XG` horizontal position (also commonly called the vorticity point).\n", "\n", "We plot the vertical integral of this quantity, i.e. the barotropic vorticity:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "zeta_bt = (zeta * ds.drF).sum(dim='Z')\n", "zeta_bt.plot(vmax=2e-4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A different way to calculate the barotropic vorticity is to take the curl of the vertically integrated velocity.\n", "This formulation also allows us to incorporate the $h$ factors representing partial cell thickness." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u_bt = (ds.U * ds.hFacW * ds.drF).sum(dim='Z')\n", "v_bt = (ds.V * ds.hFacS * ds.drF).sum(dim='Z')\n", "zeta_bt_alt = (-grid.diff(u_bt * ds.dxC, 'Y') + grid.diff(v_bt * ds.dyC, 'X'))/ds.rAz\n", "zeta_bt_alt.plot(vmax=2e-4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another interesting quantity is the horizontal strain, defined as\n", "\n", "$$ s = \\frac{\\partial u}{\\partial x} - \\frac{\\partial v}{\\partial y} \\ . $$\n", "\n", "On the c-grid, a finite-volume representation is given by\n", "\n", "$$ s = (\\delta_i \\Delta y_g u - \\delta_j \\Delta x_g v ) / A_c \\ . $$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "strain = (grid.diff(ds.U * ds.dyG, 'X') - grid.diff(ds.V * ds.dxG, 'Y')) / ds.rA\n", "strain[0,0].plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Barotropic Transport Streamfunction\n", "\n", "We can use the barotropic velocity to calcuate the barotropic transport streamfunction, defined via\n", "\n", "$$ u_{bt} = - \\frac{\\partial \\Psi}{\\partial y} \\ , \\ \\ v_{bt} = \\frac{\\partial \\Psi}{\\partial x} \\ .$$\n", "\n", "We calculate this by integrating $u_{bt}$ along the Y axis using the grid object's `cumsum` method:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "psi = grid.cumsum(-u_bt * ds.dyG, 'Y', boundary='fill')\n", "psi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that xgcm automatically shifted the Y-axis position from center (YC) to left (YG) during the cumsum operation.\n", "\n", "We convert to sverdrups and plot with a contour plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "(psi[0] / 1e6).plot.contourf(levels=np.arange(-160, 40, 5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tracer Budget Example" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "adv_flux_div = (grid.diff(ds.ADVx_TH, 'X') +\n", " grid.diff(ds.ADVy_TH, 'Y') +\n", " grid.diff(ds.ADVr_TH, 'Z', boundary='fill'))\n", "adv_flux_div" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diff_flux_div = (grid.diff(ds.DFxE_TH, 'X') +\n", " grid.diff(ds.DFyE_TH, 'Y') +\n", " grid.diff(ds.DFrE_TH + ds.DFrI_TH, 'Z', boundary='fill'))\n", "diff_flux_div" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diff_flux_div.sum(dim=['XC', 'YC']).plot.line(y='Z', marker='.', label='diffusion')\n", "adv_flux_div.sum(dim=['XC', 'YC']).plot.line(y='Z', marker='.', label='advection')\n", "plt.grid()\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Py3 pangeo", "language": "python", "name": "pangeo" }, "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }