{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# NEMO Example \n", "\n", "For this example, the NEMO output files have already been saved as netCDF with the right coordinate names. The [xorca](https://github.com/willirath/xorca) package is designed to open NEMO datasets so they are understandable by xgcm. The [xnemogcm](https://github.com/rcaneill/xnemogcm) does a similar work on idealized configurations.\n", "\n", "Below are some example of how to make calculations using xgcm.\n", "\n", "First we import xarray and xgcm:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import numpy as np\n", "import xgcm\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "plt.rcParams['figure.figsize'] = (6,10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we open the NEMO example dataset, from the BASIN configuration. The *ds* dataset contains the variables data (e.g. temperature, velocities) while the *domcfg* dataset contains the grid data (position, scale factor, etc).\n", "\n", "We get a dataset produced for the purpose of this documentation \n", "[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3736803.svg)](https://doi.org/10.5281/zenodo.3736803).\n", "In this dataset, we have not outputed the vertical scale factors in the *file_def_nemo-oce.xml* to save space. In a real case, the scale factors need to be saved to use properly the \"thickness weighted\" variables. The *domcfg* dataset has also been filtered to remove many unused variables in this notebook (e.g. *umask*)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# download the data\n", "import urllib.request\n", "import shutil\n", "\n", "url = 'https://zenodo.org/record/3736803/files/'\n", "file_name_domcfg = 'xnemogcm.domcfg.nc'\n", "with urllib.request.urlopen(url + file_name_domcfg) as response, open(file_name_domcfg, 'wb') as out_file:\n", " shutil.copyfileobj(response, out_file)\n", "file_name_basin = 'xnemogcm.nemo.nc'\n", "with urllib.request.urlopen(url + file_name_basin) as response, open(file_name_basin, 'wb') as out_file:\n", " shutil.copyfileobj(response, out_file)\n", " \n", "# open the data\n", "domcfg = xr.open_dataset(file_name_domcfg)\n", "ds = xr.open_dataset(file_name_basin)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add the vertical scale factors to *ds* (output of the model in real cases):" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "points = ['t', 'u', 'v', 'f', 'w']\n", "for point in points:\n", " ds[f\"e3{point}\"] = domcfg[f\"e3{point}_0\"]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " domcfg \n", " \n", "Dimensions: (x_c: 20, x_f: 20, y_c: 40, y_f: 40, z_c: 36, z_f: 36)\n", "Coordinates:\n", " * x_c (x_c) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19\n", " * y_c (y_c) int64 0 1 2 3 4 5 6 7 8 9 ... 30 31 32 33 34 35 36 37 38 39\n", " * x_f (x_f) float64 0.5 1.5 2.5 3.5 4.5 ... 15.5 16.5 17.5 18.5 19.5\n", " * y_f (y_f) float64 0.5 1.5 2.5 3.5 4.5 ... 35.5 36.5 37.5 38.5 39.5\n", " * z_c (z_c) float32 5.0 15.0 25.0 35.0 ... 3.38e+03 3.787e+03 4.22e+03\n", " * z_f (z_f) float32 4.5 14.5 24.5 ... 3.379e+03 3.786e+03 4.219e+03\n", "Data variables:\n", " glamt (y_c, x_c) float64 ...\n", " glamf (y_f, x_f) float64 ...\n", " gphit (y_c, x_c) float64 ...\n", " gphif (y_f, x_f) float64 ...\n", " e1t (y_c, x_c) float64 ...\n", " e1u (y_c, x_f) float64 ...\n", " e1v (y_f, x_c) float64 ...\n", " e1f (y_f, x_f) float64 ...\n", " e2t (y_c, x_c) float64 ...\n", " e2u (y_c, x_f) float64 ...\n", " e2v (y_f, x_c) float64 ...\n", " e2f (y_f, x_f) float64 ...\n", " e3t_0 (z_c, y_c, x_c) float64 ...\n", " e3u_0 (z_c, y_c, x_f) float64 ...\n", " e3v_0 (z_c, y_f, x_c) float64 ...\n", " e3f_0 (z_c, y_f, x_f) float64 ...\n", " e3w_0 (z_f, y_c, x_c) float64 ...\n", " ht_0 (y_c, x_c) float64 ...\n", " fmaskutil (y_f, x_f) float64 ...\n", "Attributes:\n", " DOMAIN_number_total: 1\n", " DOMAIN_number: 1\n", " DOMAIN_dimensions_ids: [1 2]\n", " DOMAIN_size_global: [20 40]\n", " DOMAIN_size_local: [20 40]\n", " DOMAIN_position_first: [1 1]\n", " DOMAIN_position_last: [20 40]\n", " DOMAIN_halo_size_start: [0 0]\n", " DOMAIN_halo_size_end: [0 0]\n", " DOMAIN_type: BOX\n", " nn_cfg: 2\n", " cn_cfg: BASIN\n", "\n", " ds \n", " \n", "Dimensions: (axis_nbounds: 2, t: 1, x_c: 20, x_f: 20, y_c: 40, y_f: 40, z_c: 36, z_f: 36)\n", "Coordinates:\n", " * z_f (z_f) float32 4.5 14.5 24.5 ... 3.379e+03 3.786e+03 4.219e+03\n", " * t (t) object 1050-07-01 00:00:00\n", " * x_c (x_c) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19\n", " * y_c (y_c) int64 0 1 2 3 4 5 6 7 8 ... 31 32 33 34 35 36 37 38 39\n", " * z_c (z_c) float32 5.0 15.0 25.0 ... 3.38e+03 3.787e+03 4.22e+03\n", " * x_f (x_f) float64 0.5 1.5 2.5 3.5 4.5 ... 16.5 17.5 18.5 19.5\n", " * y_f (y_f) float64 0.5 1.5 2.5 3.5 4.5 ... 36.5 37.5 38.5 39.5\n", "Dimensions without coordinates: axis_nbounds\n", "Data variables:\n", " depthw_bounds (z_f, axis_nbounds) float32 ...\n", " t_bounds (t, axis_nbounds) object ...\n", " woce (t, z_f, y_c, x_c) float64 ...\n", " deptht_bounds (z_c, axis_nbounds) float32 ...\n", " thetao (t, z_c, y_c, x_c) float64 ...\n", " so (t, z_c, y_c, x_c) float64 ...\n", " tos (t, y_c, x_c) float64 ...\n", " zos (t, y_c, x_c) float64 ...\n", " depthu_bounds (z_c, axis_nbounds) float32 ...\n", " uos (t, y_c, x_f) float64 ...\n", " uo (t, z_c, y_c, x_f) float64 ...\n", " depthv_bounds (z_c, axis_nbounds) float32 ...\n", " vos (t, y_f, x_c) float64 ...\n", " vo (t, z_c, y_f, x_c) float64 ...\n", " e3t (z_c, y_c, x_c) float64 ...\n", " e3u (z_c, y_c, x_f) float64 ...\n", " e3v (z_c, y_f, x_c) float64 ...\n", " e3f (z_c, y_f, x_f) float64 ...\n", " e3w (z_f, y_c, x_c) float64 ...\n", "Attributes:\n", " name: NEMO dataset \n", " description: Ocean grid variables, set on the proper positions\n", " title: Ocean grid variables\n" ] } ], "source": [ "print('\\n domcfg \\n', domcfg)\n", "print('\\n ds \\n', ds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the coordinates, the *_c* suffix means *center* while the *_f* suffix means *face*. Thus the (x_c, y_c, z_c) point is the T point, the (x_c, y_f, z_c) is the V point, etc." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geometry of the basin\n", "The geometry of the simulation is a closed basin, with a bottom bathymetry, going from 2000 m at the coasts, to 4000 m in the interior of the basin. Terrain following coordinates are used and the free surface is linear (fixed vertical levels).\n", "A 2 degrees Mercator grid is used." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAJDCAYAAAAFPd9oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkj0lEQVR4nO3df6xcZ33n8fcnzi8XEkg2P3Dt0FDkNoVoCWBls8tuxTZN8dKqTldNZVYFdzeSKxS6IHXVJvxT2JWlaLfQwmoTybQsTkubWlA2FgVa11vURUqTmjQQnB8bF0JyiRMnQJsgtfYm+e4fc0yGy71z5947v+5z3y9pNOc8c855Ho/P+cxznzlnTqoKSVK7Tpt2AyRJ42XQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpAlJsiHJ3yT5dDd/fpKDSR7uns/rW/amJEeTPJTkLX3lb0xyX/fah5NkqXoNekmanHcDD/TN3wgcqqqtwKFuniSvAXYCrwW2A7ck2dCtcyuwG9jaPbYvValBL0kTkGQL8NPA7/QV7wD2ddP7gGv7ym+vqhNV9TXgKHBlkk3AuVV1Z/Wudr2tb51FDRX0SV6e5BNJHkzyQJJ/PuhPDknS9/lt4NeAF/rKLq6qYwDd80Vd+Wbgsb7l5rqyzd30/PKBTh+ygR8CPldVP5/kTOAHgPfS+5Pj5iQ30vuT49cHbeTM0zbWxg3nDFmlpPXsmeeeerqqLlzNNv7Vm8+ub3/rhaUXHIEj9/2/I8A/9hXtraq9AEl+BjheVV9M8uYhNrfQuHsNKB9oyaBPci7w48AvAVTVSeBkkh3Am7vF9gGfZ4mg37jhHP7FBdctVeXUnLxsy7SbIE3MmQ/OLb3QFH3uiVu+vtptfPtbL/DJP7lgFM1Z0mWvPPaPVbVtkZffBPxskrcCZwPnJvl94Mkkm6rqWDcsc7xbfg64pG/9LcDjXfmWBcoHGqZH/8PAU8D/TPI64Iv0vlD4nj85klw0YBsrYvBK4zPJ42vWP1TGrapuAm4C6Hr0/6mqfjHJfwN2ATd3z3d0qxwA/iDJB4EfpPel691V9XySZ5NcBdwFvAP470vVP0zQnw68AfiVqroryYfovhkeRpLd9L4h5uzTXjrsaiveCf/+1WetaD1Jw3vZ355Y1vKnjuf1HvgLuBnYn+R64FHgOoCqOpJkP3A/8BxwQ1U9363zTuBjwEbgs91joGGCfg6Yq6q7uvlP0Av6xf7k+B7dGNVegJedcdGSY0n9AW9oS7NpucfmqQ8GAx+q6vP0hrqpqm8CVy+y3B5gzwLlh4HLl1PnkkFfVU8keSzJj1bVQ12j7u8eC/3JsWIrDflnf2jJ6wUkjck5X1/6nhanjuf+wF/PYT9pw5518yvAx7szbr4K/Ht6p2Z+358cKzEo4A1xabYNc4ye+jD4+1efZe9+CoYK+qq6F1jo2+QF/+QY1vxx+P6QHzbgT7zy5GqaIGkEznr0zIGvnzqez/l62bufgmF79CO3WC9+fsDPcpBfuuWpaTdB+q5H5lZ1yvmqLHWcnvogePaH8j29e+gFvr378Zpa0A9jkiFvaGutW+4+PK0Phv6w12TMdNCPg4Eu9Qw6Fqb514FGr9mgN9CllVvo+FlN+J945cklx/E1Pk3+eqUhL43epVue8thao5oKendEafw8xtaeZoZuWtz5rnnFg9Nugkbk4BOXTbsJI3XqeHMsf21Y80E/KwFvKGuQce0f0/4AMfDXhqkE/ah+z2ZSIW+Ia1YNs29O4sPg0i1PGfYzbKZ69P0XSw06h35cAW+gq0WL7dej/gCwdz+7Zirop8WA13p0ar+f9vCPxm/NnXUzyt78Na940JDXujfqY2BWvjfTi9Zc0I+KAS+9yE5P29ZU0I+ip+AOLS1unMfGLP9AYevWVNCvlgEvLW0UnaGlOmXeZ2Ky1sWXsQb86vz8ufdMuwmr8oln3jDtJqxJ17ziQb+obcSaCfqVDtu0FPJrPXCnZVrvWwsfMJ6Z04Y1E/QrMeshb3C3bbn/v7P8wbCS3v2wF1H1315Q49Fs0M9KyBvmGtagfWUWPgQcylm71kTQr6Xzcg12jcP8/WoWgl9rx5oI+uWaRm/egNckndrfZj3w/Q2c2bCuTq8ch58/9x5DXlMz6X1vVoZEtTwG/SoY8JoFdja0FIN+BTywNIvcJ7WY5oJ+3H9aejBpltkJ0UKaC/px8gDSWjHOfdVx+rWnybNuRs2A11q0Vs7M0fjNfI9+OefQj6OnYchrrXMflj36RXhwqCWj7t17lezaMvM9eknS6hj0C7A3r1a5b69PBr0kNc6gn8cej6TWGPTSOmNnZv0x6CXNjJOXbZl2E5rUTNCP4hx6ezqSWtRM0Esanp2a9cWgl6TGGfQdezjS8vjjZmuHQS9JjTPopXXKv2LXD4Nekhpn0EtS4wx6SWqcQS9JjZvJoD/xypMTrc8vpSS1bCaDXpI0Oga9tI751+z6YNBLUuMMeklqnEEvaawu3fLUtJuw7hn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glacySnJ3k7iRfSnIkyfu78vcl+UaSe7vHW/vWuSnJ0SQPJXlLX/kbk9zXvfbhJFmq/tPH88+SJPU5AfxEVX0nyRnAF5J8tnvtt6rqN/sXTvIaYCfwWuAHgT9P8iNV9TxwK7Ab+CvgM8B24LMM0ESPfjX3rvQScEnjVj3f6WbP6B41YJUdwO1VdaKqvgYcBa5Msgk4t6rurKoCbgOuXar+JoJekmZdkg1J7gWOAwer6q7upXcl+XKSjyY5ryvbDDzWt/pcV7a5m55fPpBBL0mjcUGSw32P3f0vVtXzVXUFsIVe7/xyesMwrwauAI4BH+gWX2jcvQaUD+QYvaRmffv5H+ATz7xhQrX9ydNVtW2pparq75J8HtjePzaf5CPAp7vZOeCSvtW2AI935VsWKB/IHr0kjVmSC5O8vJveCPwk8GA35n7KzwFf6aYPADuTnJXkVcBW4O6qOgY8m+Sq7mybdwB3LFW/PXpJGr9NwL4kG+h1sPdX1aeT/F6SK+gNvzwC/DJAVR1Jsh+4H3gOuKE74wbgncDHgI30zrYZeMYNGPSSNHZV9WXg9QuUv33AOnuAPQuUHwYuX079Dt1IUuMMeklqnEEvrXNeNNg+g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS44b69cokjwDPAs8Dz1XVtiTnA38EXErv5zV/oaq+PZ5mSpJWajk9+n9dVVf03UHlRuBQVW0FDnXzkqQZs5qhmx3Avm56H0PciVySNHnDBn0Bf5bki303vL24u60V3fNF42igJGl1hr3D1Juq6vEkFwEHkzw4bAXdB8NugLNPe+kKmihJWo2hevRV9Xj3fBz4FHAl8OSpG9t2z8cXWXdvVW2rqm1nnrZxNK2WJA1tyaBP8pIk55yaBn6K3p3KDwC7usV2McSdyCVJkzfM0M3FwKeSnFr+D6rqc0n+Gtif5HrgUeC68TVTkrRSSwZ9VX0VeN0C5d8Erh5HoyRJo+OVsZLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CVpzJKcneTuJF9KciTJ+7vy85McTPJw93xe3zo3JTma5KEkb+krf2OS+7rXPpwkS9Vv0EvS+J0AfqKqXgdcAWxPchVwI3CoqrYCh7p5krwG2Am8FtgO3JJkQ7etW4HdwNbusX2pyg16SRqz6vlON3tG9yhgB7CvK98HXNtN7wBur6oTVfU14ChwZZJNwLlVdWdVFXBb3zqLMuglaQKSbEhyL3AcOFhVdwEXV9UxgO75om7xzcBjfavPdWWbu+n55QMZ9JI0GhckOdz32N3/YlU9X1VXAFvo9c4vH7Cthcbda0D5QKcvtYAkrVXPPHc2B5+4bEK1/cnTVbVtqaWq6u+SfJ7e2PqTSTZV1bFuWOZ4t9gccEnfaluAx7vyLQuUD2SPXpLGLMmFSV7eTW8EfhJ4EDgA7OoW2wXc0U0fAHYmOSvJq+h96Xp3N7zzbJKrurNt3tG3zqLs0UvS+G0C9nVnzpwG7K+qTye5E9if5HrgUeA6gKo6kmQ/cD/wHHBDVT3fbeudwMeAjcBnu8dABr0kjVlVfRl4/QLl3wSuXmSdPcCeBcoPA4PG97+PQzeS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxg0d9Ek2JPmbJJ/u5s9PcjDJw93zeeNrpiRppZbTo3838EDf/I3AoaraChzq5iVJM2aooE+yBfhp4Hf6incA+7rpfcC1I22ZJGkkhu3R/zbwa8ALfWUXV9UxgO75otE2TZI0CksGfZKfAY5X1RdXUkGS3UkOJzl88oV/WMkmJEmrcPoQy7wJ+NkkbwXOBs5N8vvAk0k2VdWxJJuA4wutXFV7gb0ALzvjohpRuyVJQ1qyR19VN1XVlqq6FNgJ/O+q+kXgALCrW2wXcMfYWilJWrHVnEd/M3BNkoeBa7p5SdKMGWbo5ruq6vPA57vpbwJXj75JkqRR8spYSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glacySXJLkL5I8kORIknd35e9L8o0k93aPt/atc1OSo0keSvKWvvI3Jrmve+3DSbJU/cu6Z6wkaUWeA361qu5Jcg7wxSQHu9d+q6p+s3/hJK8BdgKvBX4Q+PMkP1JVzwO3AruBvwI+A2wHPjuocnv0kjRmVXWsqu7ppp8FHgA2D1hlB3B7VZ2oqq8BR4Erk2wCzq2qO6uqgNuAa5eq36CXpAlKcinweuCuruhdSb6c5KNJzuvKNgOP9a0215Vt7qbnlw9k0EvSaFyQ5HDfY/f8BZK8FPgk8J6qeobeMMyrgSuAY8AHTi26wPZrQPlAjtFLatbJk6fzyNyFk6ru6arattiLSc6gF/Ifr6o/BqiqJ/te/wjw6W52Drikb/UtwONd+ZYFygeyRy9JY9adGfO7wANV9cG+8k19i/0c8JVu+gCwM8lZSV4FbAXurqpjwLNJruq2+Q7gjqXqt0cvSeP3JuDtwH1J7u3K3gu8LckV9IZfHgF+GaCqjiTZD9xP74ydG7ozbgDeCXwM2EjvbJuBZ9yAQS9JY1dVX2Dh8fXPDFhnD7BngfLDwOXLqd+hG0lqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcUsGfZKzk9yd5EtJjiR5f1d+fpKDSR7uns8bf3MlScs1TI/+BPATVfU64Apge5KrgBuBQ1W1FTjUzUuSZsySQV893+lmz+geBewA9nXl+4Brx9FASdLqDDVGn2RDknuB48DBqroLuLiqjgF0zxeNrZWSpBUbKuir6vmqugLYAlyZ5PJhK0iyO8nhJIdPvvAPK2ymJGmllnXWTVX9HfB5YDvwZJJNAN3z8UXW2VtV26pq25mnbVxdayVJyzbMWTcXJnl5N70R+EngQeAAsKtbbBdwx5jaKElahdOHWGYTsC/JBnofDPur6tNJ7gT2J7keeBS4boztlCSt0JJBX1VfBl6/QPk3gavH0ShJakmSS4DbgFcALwB7q+pDSc4H/gi4FHgE+IWq+na3zk3A9cDzwH+sqj/tyt8IfAzYCHwGeHdV1aD6vTJWksbvOeBXq+rHgKuAG5K8hkWuR+pe2wm8lt53ord0oyoAtwK7ga3dY/tSlRv0kjRmVXWsqu7ppp8FHgA2s/j1SDuA26vqRFV9DThK74zHTcC5VXVn14u/jSGuYTLoJWmCklxKbzh80PVIm4HH+lab68o2d9Pzywca5stYSVqTcjKc9eiZk6rugiSH++b3VtXe72lP8lLgk8B7quqZJItta6EXakD5QAa9JI3G01W1bbEXk5xBL+Q/XlV/3BU/mWRTVR2bdz3SHHBJ3+pbgMe78i0LlA/k0I0kjVl6XfffBR6oqg/2vbTY9UgHgJ1JzkryKnpfut7dDe88m+SqbpvvYIhrmOzRS9L4vQl4O3Bf97thAO8FbmaB65Gq6kiS/cD99M7YuaGqnu/Weycvnl752e4xkEEvSWNWVV9g4fF1WOR6pKraA+xZoPwwMPTvjYFDN5LUPINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMekmagCQfTXI8yVf6yt6X5BtJ7u0eb+177aYkR5M8lOQtfeVvTHJf99qHk2Spug16SZqMjwHbFyj/raq6ont8BiDJa4CdwGu7dW5JsqFb/lZgN7C1eyy0ze9h0EvSBFTVXwLfGnLxHcDtVXWiqr4GHAWuTLIJOLeq7qyqAm4Drl1qYwa9JE3Xu5J8uRvaOa8r2ww81rfMXFe2uZueXz7Q6aNqqSTNmg0n4Zyv16SquyDJ4b75vVW1d4l1bgX+C1Dd8weA/wAsNO5eA8oHMuglaTSerqpty1mhqp48NZ3kI8Cnu9k54JK+RbcAj3flWxYoH8ihG0makm7M/ZSfA06dkXMA2JnkrCSvovel691VdQx4NslV3dk27wDuWKoee/SSNAFJ/hB4M70hnjngN4A3J7mC3vDLI8AvA1TVkST7gfuB54Abqur5blPvpHcGz0bgs91jIINekiagqt62QPHvDlh+D7BngfLDwOXLqduhG0lqnEEvSY0z6CWpcUsGfZJLkvxFkgeSHEny7q78/CQHkzzcPZ+31LYkSZM3TI/+OeBXq+rHgKuAG7rfYbgROFRVW4FD3bwkacYsGfRVdayq7ummnwUeoHfJ7Q5gX7fYPob4vQVJ0uQta4w+yaXA64G7gIu7k/fpni8aeeskSas2dNAneSnwSeA9VfXMMtbbneRwksMnX/iHlbRRkrQKQwV9kjPohfzHq+qPu+InT12+2z0fX2jdqtpbVduqatuZp20cRZslScswzFk3oXf11gNV9cG+lw4Au7rpXQzxewuSpMkb5icQ3gS8Hbgvyb1d2XuBm4H9Sa4HHgWuG0sLJUmrsmTQV9UXWPg3kAGuHm1zJEmj5pWxktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekCUjy0STHk3ylr+z8JAeTPNw9n9f32k1JjiZ5KMlb+srfmOS+7rUPJ8lSdRv0kjQZHwO2zyu7EThUVVuBQ908SV4D7ARe261zS5IN3Tq3AruBrd1j/ja/j0EvSRNQVX8JfGte8Q5gXze9D7i2r/z2qjpRVV8DjgJXJtkEnFtVd1ZVAbf1rbMog16SpufiqjoG0D1f1JVvBh7rW26uK9vcTc8vH+j0kTRVkmbQhn8sXva3JyZV3QVJDvfN762qvSvc1kLj7jWgfCCDXpJG4+mq2rbMdZ5MsqmqjnXDMse78jngkr7ltgCPd+VbFigfyKEbSZqeA8CubnoXcEdf+c4kZyV5Fb0vXe/uhneeTXJVd7bNO/rWWZQ9ekmagCR/CLyZ3hDPHPAbwM3A/iTXA48C1wFU1ZEk+4H7geeAG6rq+W5T76R3Bs9G4LPdYyCDXpImoKretshLVy+y/B5gzwLlh4HLl1O3QzeS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxSwZ9ko8mOZ7kK31l5yc5mOTh7vm88TZTkrRSw/ToPwZsn1d2I3CoqrYCh7p5SdIMWjLoq+ovgW/NK94B7Oum9wHXjrZZkqRRWekY/cVVdQyge75odE2SJI3S6eOuIMluYDfA2ae9dNzVSZLmWWmP/skkmwC65+OLLVhVe6tqW1VtO/O0jSusTpK0UisN+gPArm56F3DHaJojSRq1YU6v/EPgTuBHk8wluR64GbgmycPANd28JGkGLTlGX1VvW+Slq0fcFknSGHhlrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g15a5z7xzBum3QSNmUEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0kTkuSRJPcluTfJ4a7s/CQHkzzcPZ/Xt/xNSY4meSjJW1Zar0EvSZP1r6vqiqra1s3fCByqqq3AoW6eJK8BdgKvBbYDtyTZsJIKmwj6g09ctuJ1vfxb0pTtAPZ10/uAa/vKb6+qE1X1NeAocOVKKmgi6CVpjSjgz5J8McnuruziqjoG0D1f1JVvBh7rW3euK1u201fYWEmaefnHk5z54Nykqrvg1Lh7Z29V7Z23zJuq6vEkFwEHkzw4YHtZoKxW0jCDXpJG4+m+cfcFVdXj3fPxJJ+iNxTzZJJNVXUsySbgeLf4HHBJ3+pbgMdX0jCHbiRpApK8JMk5p6aBnwK+AhwAdnWL7QLu6KYPADuTnJXkVcBW4O6V1G2PXpIm42LgU0mgl71/UFWfS/LXwP4k1wOPAtcBVNWRJPuB+4HngBuq6vmVVGzQSxqrR+YunHYTZkJVfRV43QLl3wSuXmSdPcCe1dbt0I0kNc6gl9YxryNZHwx6SWqcQS9JjZvJoD/r0TMnWp9/vkpq2UwGvSRpdAx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfTSOuVpxeuHQd9xp5eWZzW38NRkGfSS1Lhmgn4UvQt79Vov3NfXl2aCXtLaN8H7u64rBr20ztibX38M+nk8CCS1xqBfgGGvVrlvr08GvbROjDLkPbVybTHoF2HPR1Irmgr6UfcyDHu1wn15fZv5oH9k7sKp1u8BorXOfVgzH/SzwANFa5X7rsCgH5oHjNYa91mdYtBLDTLk1c+gXwYPHslTK9cig36ZDHvNOvdRzddc0E+it+GBpFnlvqmFnD7tBqxV8w+onz/3nim1ROvZpIN9uR2paZ8erZ41EfSPzF3IpVueGnr5g09cxjWveHCMLfp+ix1wfgBoFOypazXWRNCvxDTCfiGDDlA/BNRv1sPc3vzatWaCfrm9epidsF/MSg9sPyBm26wH9kqM87uvl/3tibFtWz1rJuhh5WEPzHTgL1eLQaLZtNKAtzc/W9ZU0K9Gi4Evjcu4z1475+s11u3re6250ytX21PwYg9pcQefuGzVx4i9+dmz5oJ+FEaxM0utGfcxcdajZ451+1rcmhy6WclY/UIczpFGG/D25mfTmgz6UTPwtR75V+36sWaHbsbRc3DH13owrqFLe/Oza6Z69Od8vXj2hwL0xvNOvPLkxNuw0AFgT19rmR0YrSrok2wHPgRsAH6nqm5e7jZe9rcn+PtXn7Wi+kc1Vr+U+QeKwa9ZNa1Qtzc/21Yc9Ek2AP8DuAaYA/46yYGqun+pdc98cI6Tl21ZadVTZ69fs2BWeuqG/OxbTY/+SuBoVX0VIMntwA5gyaAfpUn16peynIPODwUNMisBrnasJug3A4/1zc8B/2x1zVkfRnkg+6ExG9ZrONubXxtWE/RZoOz7rmtOshvYDXD2aS9dRXWLm5Ve/TSs14CRNLzVBP0ccEnf/Bbg8fkLVdVeYC/Ay864aGw/cDG/Z7Feg18al9X03vuviu3/nRt/uXIyVhP0fw1sTfIq4BvATuDfDbty/xey/WfejOoUy5XslH44aL2Y5JDLsD99cOaDc2Nuyfq14qCvqueSvAv4U3qnV360qo4sZxvDhv184zq/vrXxRj+4Rqu1/WNcBgX7Yr15Q368VnUefVV9BvjMarYxTNjPN2hHmsZFVrPKYNK4rOQHygz56ZmpK2PnGxT2i1nJDuiHg9azSfyqpCE/XTMR9Iv16ifFn08dbC1/EPp/O31++Tp9MxH0MHgIZ5Dl9vi1fIallrKSO0bZm5+cmQl6WFnPfrk7mB8M0tLGcas/h2ymZ6aCHsY/jOO9KqXJM+Sna+Z/j94xPWltM+Snb+Z69PD9v2653LCf9Je50nphx2ttmsmgh9X9lLE7ozR77M1Pz0wP3bhjSG3wWJ6umQ566O0g7iTS2uTx+6Ik25M8lORokhsnWffMDt3M584iaa1azR35RmHme/SS1IDv3pGvqk4Cp+7INxEGvSSN30J35Ns8qconOnTzzHNPPf25J275+iIvXwA8Pcn2DGBbFjdL7bEtC5ultsDK2/NDq634meee+tPPPXHLBavdzpDOTnK4b35vd+MlGPKOfOMy0aCvqkV/NzfJ4araNsn2LMa2LG6W2mNbFjZLbYHptqeqtk+j3gUMdUe+cXHoRpLG77t35EtyJr078h2YVOVr5qwbSVqrRnFHvtWYpaDfu/QiE2NbFjdL7bEtC5ultsDstWcqRnFHvpVKlb/mKEktc4xekho39aCf5mXBi7TnkST3Jbl33qlSk6j7o0mOJ/lKX9n5SQ4mebh7Pm+KbXlfkm907829Sd46obZckuQvkjyQ5EiSd3fl03pvFmvPxN+fJGcnuTvJl7q2vL8rn/h7M6AtU9lv9KKpDt10lwX/X/ouCwbeNqnLghdp0yPAtqqa+HnISX4c+A5wW1Vd3pX9V+BbVXVz90F4XlX9+pTa8j7gO1X1m+Ouf15bNgGbquqeJOcAXwSuBX6J6bw3i7XnF5jw+5MkwEuq6jtJzgC+ALwb+LdM+L0Z0JbtTGG/0Yum3aOf6mXBs6aq/hL41rziHcC+bnofvUCZVlumoqqOVdU93fSzwAP0riqc1nuzWHsmrnq+082e0T2KKbw3A9qiKZt20E/1suBFFPBnSb6YZPeU2wJwcVUdg17AABdNuT3vSvLlbmhnIkMl/ZJcCrweuIsZeG/mtQem8P4k2ZDkXuA4cLCqpvbeLNIWmPJ+s95NO+inelnwIt5UVW8A/g1wQzeEoZ5bgVcDVwDHgA9MsvIkLwU+Cbynqp6ZZN1Dtmcq709VPV9VV9C72vLKJJdPot5ltGWq+42mH/RTvSx4IVX1ePd8HPgUveGlaXqyGxM+NTZ8fFoNqaonuwP5BeAjTPC96cZ8Pwl8vKr+uCue2nuzUHum+f509f8d8Hl6Y+JT3W/62zLt90XTD/qpXhY8X5KXdF+ukeQlwE8BXxm81tgdAHZ107uAO6bVkFPB0fk5JvTedF/y/S7wQFV9sO+lqbw3i7VnGu9PkguTvLyb3gj8JPAgU3hvFmvLtPYbvWjqF0x1p1r9Ni9eFrxnim35YXq9eOhdNfwHk2xPkj8E3kzv1/6eBH4D+F/AfuCVwKPAdVU19i9JF2nLm+n9+V3AI8AvnxoHHnNb/iXwf4D7gBe64vfSGxefxnuzWHvexoTfnyT/lN6XrRvoddz2V9V/TvJPmPB7M6Atv8cU9hu9aOpBL0kar2kP3UiSxsygl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcf8f1RZvavMl+oEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.contourf(domcfg.glamt, domcfg.gphit, domcfg.ht_0)\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating the grid object\n", "\n", "Next we create a `Grid` object from the dataset.\n", "All the axes are here non-periodic.\n", "The `metrics` dict contains the scale factors." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Z Axis (not periodic, boundary=None):\n", " * center z_c --> left\n", " * left z_f --> center\n", "X Axis (not periodic, boundary=None):\n", " * center x_c --> right\n", " * right x_f --> center\n", "Y Axis (not periodic, boundary=None):\n", " * center y_c --> right\n", " * right y_f --> center\n" ] } ], "source": [ "\n", "metrics = {\n", " ('X',): ['e1t', 'e1u', 'e1v', 'e1f'], # X distances\n", " ('Y',): ['e2t', 'e2u', 'e2v', 'e2f'], # Y distances\n", " ('Z',): ['e3t_0', 'e3u_0', 'e3v_0', 'e3f_0', 'e3w_0'], # Z distances\n", " #('X', 'Y'): [] # Areas TODO\n", "}\n", "grid = xgcm.Grid(domcfg, metrics=metrics, periodic=False)\n", "print(grid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that xgcm identified three different axes: X (longitude), Y (latitude), Z (depth)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computation examples\n", "### Horizontal gradient of SST\n", "\n", "We will compute the horizontal component of the gradient of SST in the longitude\n", "direction as a first example to understand the logic behind the\n", "xgcm grids.\n", "\n", "We want to compute\n", "$\\frac{\\partial SST}{\\partial x}$.\n", "The SST is the variable `tos` (Temperature Ocean Surface) in our dataset.\n", "\n", "In discrete form, using the NEMO notation, the derivative becomes [1]\n", "\n", "$$\\frac{\\partial SST}{\\partial x} = \\frac{1}{e_{1u}} \\delta_{i+1/2} SST\n", "= \\frac{1}{e_{1u}} (SST_{i+1} - SST_i) \\ .$$\n", "\n", "The last T point is an earth point here, such as the 2 last U points: we set up the\n", "`boundary` argument to 'fill' and the fill value to zero (this value does not play an important role here, as we fill earth points).\n", "\n", "The gradient is first computed with the `diff` function and\n", "then with the `derivative` function, the result is the same as the `derivative`\n", "function is aware of which scale factor to use.\n", "\n", "
\n", "[1] NEMO book v4.0.1, pp 22" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Coordinates:\n", " * t (t) object 1050-07-01 00:00:00\n", " * y_c (y_c) int64 0 1 2 3 4 5 6 7 8 9 ... 30 31 32 33 34 35 36 37 38 39\n", " * x_f (x_f) float64 0.5 1.5 2.5 3.5 4.5 5.5 ... 15.5 16.5 17.5 18.5 19.5\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray ()>\n",
       "array(True)
" ], "text/plain": [ "\n", "array(True)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "grad_T_lon0 = grid.diff(ds.tos, axis='X', boundary='fill', fill_value=0) / domcfg.e1u\n", "grad_T_lon1 = grid.derivative(ds.tos, axis='X', boundary='fill', fill_value=0)\n", "print(grad_T_lon1.coords)\n", "(grad_T_lon0 == grad_T_lon1).all()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected the result of the 2 operations is the same.\n", "The position of the derivative is now on the U point." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Divergence Calculation\n", "\n", "Here we show how to calculate the divergence of the flow.\n", "The precise details of how to do this calculation are highly model- and configuration-dependent (e.g. free-surface vs. rigid lid, etc.)\n", "In this example, the flow is incompressible, without precipitations\n", "or evaporation, with a linear free surface, satisfying the continuity equation\n", "\n", "$$ \\vec{\\nabla} \\cdot \\vec{u} = \n", " \\frac{\\partial u}{\\partial x} + \\frac{\\partial v}{\\partial y}\n", " + \\frac{\\partial w}{\\partial z} = 0 \\ .$$\n", " \n", "In non-linear free surface, the divergence of $\\vec{u}$ is not 0 anymore, it is linked to\n", "the time variation of the $e_{3}$ scale factors.\n", "\n", "In discrete form, using NEMO notation, the equation becomes [1]\n", "\n", "$$ \\vec{\\nabla} \\cdot \\vec{u} = \\frac{1}{e_{1t}e_{2t}e_{3t}} \\left[\n", " \\delta_i(u \\cdot e_{2u} \\cdot e_{3u})\n", " + \\delta_j(v \\cdot e_{1v} \\cdot e_{3v}) \\right]\n", " + \\frac{1}{e_{3t}} \\delta_k(w) \\ .$$\n", "\n", "The equation for the divergence calculation could be simplified due to the geometry of our basin, but we will keep it in the general form.\n", "\n", "
\n", "[1] NEMO book v4.0.1, pp 22" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3d flow\n", "\n", "We can't use the `grid.derivative` function for *w* even if this is a simple vertical derivative: xgcm does not use the temporaly varying scale factors but the *domcfg.e3t_0* initial vertical scale factors. In this simple linear free surface case, the scale factors don't vary in time, but we will use here a more general approach.\n", "\n", "For *u* and *v* we use the `grid.diff` as this is not a simple derivative, even if the horizontal scale factors are fixed in time in any case.\n", "The first T point along x is an earth point, we can use a 'fill' `boundary` with a 0 `fill_value`. The same argument applies along y and z.\n", "\n", "⚠ We need to use a negative sign for the vertical derivative, as the k-axis increases with depth.\n", "\n", "⚠ If the model is ran without linear free surface the e3 scale factors are not constant in time, we need to take the e3 from *ds* and not from *domcfg*. The data are are \"thickness weighted\". To stay general, we use here the e3 scale factors from *ds*." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray ()>\n",
       "array(5.44218669e-19)
" ], "text/plain": [ "\n", "array(5.44218669e-19)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bd={'boundary':'fill', 'fill_value':0}\n", "\n", "div_uv = grid.diff(ds.uo * domcfg.e2u * ds.e3u, 'X', **bd) / (domcfg.e1t * domcfg.e2t * ds.e3t) \\\n", " + grid.diff(ds.vo * domcfg.e1v * ds.e3v, 'Y', **bd) / (domcfg.e1t * domcfg.e2t * ds.e3t)\n", "\n", "div_w = - grid.diff(ds.woce, 'Z', **bd) / ds.e3t\n", "\n", "div_uvw = div_uv + div_w\n", "div_uvw.max()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected the divergence of the flow is zero (if we neglect the truncation error)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Vertical velocity\n", "\n", "In NEMO the vertical velocity is computed from the divergence of the horizontal velocity.\n", "In non-linear free surface, the vertical velocity can't be computed offline because it also takes the time\n", "variations of the vertical scale factors into account.\n", "However, we are using here a linear free surface, so that\n", "\n", "$$\n", "w(z) = \\int_{bottom}^z \\vec{\\nabla}_h \\cdot \\vec{u} \\, \\text{d}z' \n", "= \\int_{surf}^z \\vec{\\nabla}_h \\cdot \\vec{u} \\, \\text{d}z' - \\int_{surf}^{bottom} \\vec{\\nabla}_h \\cdot \\vec{u} \\, \\text{d}z' \\ .\n", "$$\n", "\n", "This is written in discrete form\n", "\n", "$$\n", "w(n) = \\sum_0^n \\left(\\vec{\\nabla}_h \\cdot \\vec{u}(k)\\right) e_{3t}(k)\n", "- \\sum_0^{n_{bot}} \\left(\\vec{\\nabla}_h \\cdot \\vec{u}(k)\\right) e_{3t}(k) \\ .\n", "$$\n", "\n", "We use the `grid.cumsum` to perform the integration, and then we remove the total integral.\n", "This is shown here with 2 methods that give the same numerical result, however the first method runs faster." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray ()>\n",
       "array(True)
" ], "text/plain": [ "\n", "array(True)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w = grid.cumsum((div_uv*ds.e3t), axis='Z', boundary='fill', fill_value=0) # integral from top\n", "w = w - w.isel({'z_f':-1}) # now from bot\n", "\n", "w_alt = grid.cumsum((div_uv*ds.e3t), axis='Z', boundary='fill', fill_value=0) - (div_uv*ds.e3t).sum(dim='z_c')\n", "\n", "(w == w_alt).all()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remark: the first method runs faster, indeed we perfom only once the integral." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.19 ms ± 86.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n", "6.41 ms ± 232 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "def w1(grid, div_uv, ds):\n", " w = grid.cumsum((div_uv*ds.e3t), axis='Z', boundary='fill', fill_value=0) # integral from top\n", " return w - w.isel({'z_f':-1}) # now from bot\n", "\n", "def w2(grid, div_uv, ds):\n", " return grid.cumsum((div_uv*ds.e3t), axis='Z', boundary='fill', fill_value=0) - (div_uv*ds.e3t).sum(dim='z_c')\n", "\n", "%timeit w1(grid, div_uv, ds)\n", "%timeit w2(grid, div_uv, ds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now compare the computed vertical velocity with the one outputed by the model, at the bottom of the uper grid cell." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAJcCAYAAAD9zQMxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABad0lEQVR4nO3de5xkVX3v/c+3q7vnDgzMCAOMoEiiQmQgc8D7QUUDhIga9YiJwcsRSSTqMZ7IgyfCSUyCxHs06KgjYBDv6DweQAkR8QYPIwKCoCIHdGBkGO7DXHq66/f8sXfPVNdUde+9uu7zfb9e+9VVtffaa+1d1bVqr7X2bykiMDMzmzTU7QKYmVlvccVgZmZTuGIwM7MpXDGYmdkUrhjMzGwKVwxmZjaFKwYzM5vCFcMsSJon6f+V9Iikr3S7PGZmreCKYXZeCewL7BMRr2rVTiW9QNJ38wrnrmm2+6+SQtL7ptnmHEnbJW2qWZ5csByvl/SDhEMosu8zJK2VtE3SBQ3Wv0jS7ZI25+fioJp10x6TpIPzNJvzfRw3Q1mabi/prLp8tkiqSloyzf5eK+luSY9L+oakvWvWvVrSj/K8ri5wnqY7D5L0fkkP5Mt5kpRynDOVu8G+5khaLelRSb+T9M669Ssk/STP6yeSVsx0rNY7XDHMzkHALyNivMX7fRxYDfzPZhtIGgE+ClxXYH9fioiFNcudLSrnbNwLvI/sOKfIv3S/DvwdsDewFvhS3WbTHdMlwE+BfYD3AF+VtHSasjTdPiL+qTYf4P3A1RGxsdGOJB0GfAp4HdmPhs3Av9Vs8iDwEeDcacozua+ZzsNpwMuAI4BnACcBb0k5zgLlrncOcCjZ/8ALgL+VdHy+r1Hgm8C/A4uBC4Fv5q9bP4gIL1lYkP8GbKpZtpF9ATTb/n8DY8D2fPs3taFMxwF3NVl3JnAecAHwvmn2cQ7w7wl5Pw3YCkzkx/dwm877+4AL6l47DfhRzfMFwBbgqTMdE/B7+Xu3qOa17wOnz3Z7QMCvgVOnOZ5/Ar5Q8/yQ/HOyqG67/z7d56vgefgRcFrN+jcB16YcZ9Fy16y/B3hJzfN/AL6YP35Jvl41638DHN+Oz5CX1i++YshFxI5foMD+wJ1kv7CabX822T/TZLrP1m+TX5o/PM3yxJSy5s0JbwT+vmCSP5H0oKRbJf1lkQQRcRtwOvDj/Pj2alKWf5vm+G4uWL56hwE31ZTlcbIv5MMKHNNhwJ0R8VjNazfVpa3Pq+j2zyP7Nf21EmX/NdkX7O9Nk6bovurPw5T19eWW9C1JZ9ZsO91xTltuSWdK+lb+eDHZ/0izvA8Dbo68RsjdTPP3oG3y5q4Nkm5p0f4mJN2YL2tasc9eNNztAvQaSUPAF8h+zX1qNvuKiC/k+2q1jwF/FxGbpmlSnvRlYBVwH3AM8DVJD0dE00qvjIj4K+CvWrGvGguB++teewRYlD+e7pgW5tvWpz1gmryKbn8q8NWI2DRD2Rvtb1GDbWcy03moz+sRYKEkReakAuU6YIb1iwAiorbpa2HN+iLlql/fSRcAHwcuatH+tkTEihbtq2f5imFX/0j2AX5btwvSiKQ/Ibu8r29zbygifh4R90bERET8iKxf4pVtLeTsbQL2qHttD+AxmPGYpk2bX2FMdiQ/b6btJ0maB7yKrL188rXn1ezr1iL5lzTTvurX7wFsqvulnrqv+vX1+5pcP9t9tVVEXEPWr7ODpEMkXZF3in9f0lM7Xa5e54qhhqTXAKcAr4yI7S3Y35/VjWipX1Kakl4ErMxHgvyOrG/kHZK+WTB9kLWVF912WpI+Oc3x3TpT+iZuJetQncxjAVmbd7P91R7TrcCTJdX+Oj1iMm1EHBY7O5O/P9P2NV5B9gVz9Y5MI75fs6/JZpL6sj8ZmAP8csaj3tVM52HK+iblrt3XdMdZuNwR8RCwfpq8bwWeUTdC6hnTlK3TVgF/HRF/CLyL6TvZ681VNpruWkkva0vpekG3Ozl6ZQGOJLtsX1EizTkkdOwW2O8QMBc4Abg7fzyar1sE7FezfAn4MLB3k32dTDYyRMDRZJ2Cp9asvxo4p0na44G7JvNu8TEO58f1z8Dn88fD+bqlZE0Pf5q//n5qOlULHNO1wAfytC8HHgaWTlOWGbcHvgP8fYHjOgx4lKw/YgHZyJwv1qyv5PmcDlyTPx5psq+ZzsPpwG1kzUH7k33xNuxkn+k4Zyp3g32dC3wvfx+eSlZRHJ+vG80/t28nq1zOyJ+3/HNU8LN2MHBL/nghWQf+jTXLbfm6VwC3NFi+XbOv/fO/T87/Nw7pxjG1/Zx1uwC9spB9yY8zdWTS5QXStKNiOJbsV3DtcnWTbS+gZlRS/o+9qeb5JcAD+fHcDrytLv2vgRc32fco8H/IfilvbMP5rj/Gc2rWH5eXdwtZ5XVwiWM6OE+zBfgFcNwMZZl2e7Iv3nHgKQWP7bVko3AeJxu2uXfNutc3OO4LptnXdOdBZCPTHsyX85g6Euhy4KwSxzlduc+q/X8g+8JfTVaZ3Ae8s25fRwI/yfO6ATiy1f8nJT5rB7OzYtgDWN+i/V5A1rrQleNq56L8AG03JOlA4CsR8axul8WsXSQdDHwrIg7Pn/8I+HBEfCVv7npGRNw03T7ydIuBzRGxLb/H5MfAyRHx8zYWvytcMZjZwJJ0CdkV+BKyK5uzgf8EzgeWASNkTWYzDv2W9GyymwCrZM29H4kGw9QHgSuGGeQdqAc1WPWWiLi40+UxM2s3VwxmZjZFX9zgtmTJkjjoiUk3CZvZbuaGn/50Y0RMFxtrRss1L7ZSLbz9Rsa+HRHHzybPXtIXFcNBT3wiP/zhD7tdDDPrA/Pmz797tvvYSpU/ZVnh7T/F3U2j7fajvqgYzMw6SUCl6G2gUOBW0P7iisHMrE5WMZSoGVwxmJkNvlJXDAPGFYOZWZ3SVwwDxhWDmVkdCUaHXDEMnMqjvyudpjp/cVpmUXxYW99QYuDdlHORklfqOe9kXikSyqdq2syyURkpn9fWtMjZY/P3KZ1muIvfy6U7nwfMwFYMZmbp5KYkMzPbSezek9W4YjAza8BXDGZmtoPkPgYzM6vjK4YBVF1YPnTJdqWdjpTRE9XC0y7vNJR4e2Xsxh/wfpMS7Dj1c7F1ony6eaPzk/LavL38yK49RrvXyu9RSWZmNoVvcDMzs134isHMzHbIOp9335rBFYOZWQO+YjAzsx3c+TyoqhOlk4ykhgdKuEeykhDfJobS3q6h6gDGcuqgSIwblRLDKCqjpdMMjT1eOg3AnDkLS6fR1rS4TPNG5yWl6xY5JIaZmdXzFYOZme3gzmczM5vCfQxmZrYLXzGYmdkOvmJoE0lzgWuAOXk+X42IsyWdA7wZuD/f9KyIuKwNBSidpJowoxWAEgLcpIw+SZU6qqaXpcZ/SnqvOvjLcTwh7NHI8JykvLaNlx+tNj/xs7QlIVbSSBdjJYGvGNplG/DCiNgkaQT4gaTL83UfjogPtDFvM7NkEgy5Ymi9iAhgU/50JF/SwkCamXWU0G7cltTWazVJFUk3AhuAKyPiunzVGZJulrRa0uImaU+TtFbS2vs3bmxnMc3MphIMVVR4GTRtrRgiYiIiVgAHAkdLOhw4HzgEWAGsBz7YJO2qiFgZESuXLik/t4KZWSoBqgwVXgZNR44oIh4GrgaOj4j78gqjCnwaOLoTZTAzK0ygigovg6ZtFYOkpZL2yh/PA44Dbpe0rGazlwO3tKsMZmZJVLwZaRCbkto5KmkZcKGkClkF9OWI+Jakz0taQdYRfRfwlrbknjAsUQmB96w7lPibJiWwHYnBC1PyqqQMPR3bXj4NMDxc/rgiEs9F9N+Xp4Za97tZ0mrgJGBDRBzeYL2AjwInApuB10fEDS0rQEntHJV0M3Bkg9df1648zcxaQXnncwtdAHwcuKjJ+hOAQ/PlGLK+2GNaWYAyfOezmVk9iaHRSst2FxHXSDp4mk1OBi7Kh/lfK2kvScsiYn3LClGCKwYzszqi9BXDEklra56viohVJdIfAPy25vm6/DVXDGZmPUGgoVIVw8aIWDm7HHfRtRuCXTGYme1CDHX2/oR1wPKa5wcC93ayALVcMdSIobQ2xZTAbEkicYpOB9HbKeFcJAfsS8kr5aOUGPwxSeJnMJIOrIsjmfL7GDpoDVlEiC+SdTo/0q3+BXDFYGa2C7W4YpB0CXAsWV/EOuBssvhxRMQngcvIhqreQTZc9Q0tyzyBKwYzswZa2ZQUEafMsD6At7Ysw1lyxWBmVk+DGeqiKFcMZmZ1BAyVG5U0UFwxmJnVEwMZNbUoVww1HCupRidHQCXklRwrKSGvoHV3wM4kaQBU8nuVlixFYgm7ahCD4xXlisHMrJ77GMzMrJbclGRmZvXclGRmZjuVj5U0UFwxmJnVUedjJfWUwa0YUkZq9HpModTRJyl6/Vwkih4/rrSQQr19TH2p87GSesrgVgxmZqnc+WxmZlOppXM+9xtXDGZmdbI5n10xmJnZDnJTkpmZ1ejTPgZJRxXYbHtE/Gy6DdpWMUiaC1wDzMnz+WpEnC1pb+BLwMHAXcCrI+KhluefFBMnkUdA7ZRyXJ08Fx3MK2UEVOrEdL2uU5Mctk7f9jF8D7ie6SNhPYns+7epdl4xbANeGBGbJI0AP5B0OfAK4KqIOFfSmcCZwLvbWA4zs3IEqnQueGILXR8RL5xuA0n/OdNO2lYlRmZT/nQkXwI4Gbgwf/1C4GXtKoOZWQrlfQxFl14xU6VQdJu2HpGkiqQbgQ3AlRFxHbDv5CTX+d8nNEl7mqS1ktbev3FjO4tpZjaVYGhoqPDSayQ9R9KC/PGfS/qQpIOKpm/rEUXERESsAA4EjpZ0eIm0qyJiZUSsXLpkSdvKaGbWSD9eMdQ4H9gs6Qjgb4G7gYuKJu7IEUXEw8DVwPHAfZKWAeR/N3SiDGZmhak/m5JqjEfEZNP9RyPio8CioonbOSppKdmwqIclzQOOA94PrAFOBc7N/36zXWXomE6NdBnQ0Tu9TqlDajp0DlPjPw2lDIFKzKvaZ6OSJDE00tej+R+T9P8Afw48X1KFrJ+3kHYe+TLgwrxAQ8CXI+Jbkn4MfFnSm4DfAK9qYxnMzMrr0/sYavw34LXAmyLid5KeCPxL0cRtqxgi4mbgyAavPwC8qF35mpm1Qj+GxJD0beAK4PKI+NDk6xHxG0r0MfT1tZKZWTtIfXuD26lkfbnnSPo94DqyiuKqmtsHZuSKwcysgX5sSoqI3wEXABdIGgKOAU4A/lbSFuA7EXHeTPtxxWBmVk/9H0QvIqrAj/PlvZKWAH9UJK0rBjOzBvq0KQkASU8C3gYcRM33fES8tEj6ga0YYqj8oUViBLOU0YxD6SH7Bs8gBiEk7fOUPDTWWkoSQ/0ZK2nSN4DPkt0eUPofbGArBjOz2ejzpqStEfGx1MSuGMzM6vX/fQwflXQ28B2ySNcARMQNRRK7YjAz20XfDled9AfA64AXsrMpKfLnM3LFYGZWR/1/xfBy4MkRMZaS2BWDmVm9/h+uehOwF4lBSge2YojhOaXTpAdLSxjNVO3xkTidnNozQeoIMlUnyuc1lDY6pWMjjBLPedIpHNS5Rxvo86akfYHbJV3P1D6G3Xu4qplZMgkl/iBovksdD3wUqACfiYhz69YfSxZt+v/mL309Iv4+MbuzE9MBrhjMzBprYcWQR5n+BPBiYB1wvaQ1EfHzuk2/HxEnzTa/iPjebNL39bWSmVl7CIaGii8zOxq4IyLuzDuEv0g2iU5rSy19qxXb+IrBzKyeQOXufF4iaW3N81URsarm+QHAb2ueryMLcFfvWZJuAu4F3hURt5YpBPBcSWumWS/g6TPtxBWDmdkuVLYpaWNErJx+h7uoH51wA3BQRGySdCJZWItDyxSCYlchMw5hHdiKYetE+REhc4fSRpEE5dsiU8Z2dHIkjqcRbYEejwHVybBMQ/02mEm0tI+B7Aphec3zA8muCnaIiEdrHl8m6d8kLYmIjUUzmW3fwqQB/Y80M0un/M7noksB1wOHSnqSpFHgNWQB7nbmKe0nZb/+JB1N9v38QIsPrZCBvWIwM0vW4iuGiBiXdAbwbbLhqqsj4lZJp+frPwm8EvhLSePAFuA1Ed0Jt+uKwcxsF6X7GGYUEZcBl9W99smaxx8HPt6KvCSdBFyWT9ZTmpuSzMwaaHFTUqe9BviVpPMkPa1sYl8xmJnVU+uvGDopIv5c0h7AKcDnJAXwOeCSiHhspvRtqxgkLQcuAvYjC/u6KiI+Kukc4M3A/fmmZ+WXWC01f3xT6TTV0QVJeVUTmgErE9vL56O0D+pQdbx8XgmxpgCGxrfNvFGd7UOjpdMMR8JIK9LiHqXH0OrxuFEtLsdg6e+KAbJRTpK+BswD3kEWcfV/SvpYRPzrdGnbecUwDvxNRNwgaRHwE0lX5us+HBEfaGPeZmbpyt/g1lMkvRR4A3AI8Hng6IjYIGk+cBvQnYohItYD6/PHj0m6jezuPzOz3ibB8Ei3SzEbryT7AX5N7YsRsVnSG2dK3JFrXUkHA0cC1+UvnSHpZkmrJS1ukuY0SWslrb1/Y+H7O8zMWiCLrlp06UHr6ysFSe8HiIirZkrc9opB0kLga8A78jv7zie7vFlBdkXxwUbpImJVRKyMiJVLlyxpdzHNzKZqbRC9Tntxg9dOKJq4raOSJI2QVQoXR8TXASLivpr1nwZmjPRnZtZRbZiPoRMk/SXwV8Ahkm6uWbUI+GHR/bRzVJKAzwK3RcSHal5flvc/QNZLfks78q88sn7mjepsW/p7aXkljAqJSvn2y5TRTwCVhNEx49W0vMqPL9o1klgRSpxhrhopo5J6eza7ftB/E7/17aikLwCXA/8MnFnz+mMR8WDRnbTziuE5wOuAn0m6MX/tLOAUSSvIvg/uAt7SxjKYmZUnerWJaCYREXdJemv9Ckl7F60c2jkq6Qc0Hird8nsWzMxaSahfh6t+ATgJ+AnZj+/a7+AAnlxkJ77z2cysXuvDbnfE5LSgEfGk2eynL6+VzMzaK+9jKLr0GEkvl7RnzfO9JL2saHpXDGZmDfR5EL2zI+KRyScR8TBwdtHEbkoyM6vX50H0aPyjv/D3/eBWDB2sxYfHygfsi4QgdUOpH9Rq+eGWlV6fizFxCGnSsMnUqVISypgS5K+jEofgTiSO+O2q/h5uvFbSh4BPkH2C/5qsQ7qQvj5yM7P2UFYxFF16z18DY8CXgK8AW4FdhrA2M7hXDGZmsxC9+YVfSEQ8DpyZz8lQjYhSzRr9e+RmZu0i+vqKQdIfSPop8DPgVkk/kXR40fS+YjAz24X6MY5HrU8B74yI7wJIOhZYBTy7SGJXDGZmjfTmMNSiFkxWCgARcbWkwlNUumKo8fj2tKETo0PlT+N2lU9T6YcfMAmX1UmHlXj5nhKHMPW0d3Qa0QSpU4KmqPTZd2zQ330MwJ2S/o5s9jaAPwf+b9HEfX3kZmZtob4flfRGYCnwdeDS/PEbiib2FYOZWSO9+YVfSEQ8BLwtNb0rBjOzXagvKwZJ/y/T3JIZES8tsh9XDGZmDfRpH8MHWrETVwxmZo30YcUQEd+bfCxpHvDEiPhF2f0MbMUQlfKTTO41mvZB0Lax0mkmhsrHShrd+mjpNABbRhaVTjO3uj0prxQpcZkiOvfR7flpRDv4BVZNmJIWgPGJ1hak3dTf9zFI+hOyq4dR4En5rJl/X7Qpqf+qRDOzTujvUUnnAEcDDwNExI3AwUUTD+wVg5nZbPRpH8Ok8Yh4RIlXPa4YzMx2oX6/8/kWSa8FKpIOJRu6+qOiifv6yM3M2qLPg+iRhd0+DNgGfAF4BHhH0cS+YjAz24UgIdRND/n9iHgP8J6UxG07cknLgYuA/YAqsCoiPippb7LJIw4G7gJend+l19r8H1xXOs33HluclNdxc9eXTvPAgnml0yyfeLx0GoAYXlg6jSa2JeU1NlI4TtcOwx2MD9TJgSYTCcdVSRiVVE2c9S0lLpOq40l5DfVh40Sf9zF8SNIyskl6vhgRt5ZJPOORS/onSXvVPF8s6X0F9j0O/E1EPA14JvBWSU8HzgSuiohDgavy52ZmvaPPYyVFxAuAY4H7gVWSfibpfxVNX+SIToiIh2syfAg4sUDB1kfEDfnjx4DbgAOAk4EL880uBF5WtLBmZh0zeS9DkaUHRcTvIuJjwOnAjcB7i6YtUjFUJO24Gyu/m67U3VmSDgaOBK4D9o2I9ZBVHsATmqQ5TdJaSWvv37ixTHZmZrPU+isGScdL+oWkOyTt0lKizMfy9TdLOiq59NLTJJ0j6Rbg42Qjkg4smr5IH8O/A1dJ+hxZcKY3svMXf5ECLgS+BrwjIh4tOq42IlaRzTjEHx51VOcaoc3MaG0fg6QK8AngxcA64HpJayLi5zWbnQAcmi/HAOfnf1N8DrgEeElE3Fs28YwVQ0ScJ+lm4DiyQVz/EBHfLrJzSSNklcLFEfH1/OX7JC2LiPV558iGsoU2M2u71vYdHA3cERF3Akj6Ilmzem3FcDJwUUQEcK2kvSa/K8tmFhHPnE1hC41KiogrgCsarZP044h4VoPXBXwWuC0iPlSzag1wKnBu/vebZQtdRCzev3SaY/YsP3oHIB6dWzrNXnPKjySJbWlxah7fXv6Ca36Uj/8EEAmjklJGukQHhxL2+uiU1EFdKS3jE0Npn8Eq/RUrKaSyM9wtkbS25vmqvNVj0gHAb2uer2PXq4FG2xwAlB/2OEut+O9q9q34HOB1wM8k3Zi/dhZZhfBlSW8CfgO8qgVlMDNrnShd4W6MiJXTrG9Uy9TnUGSbjmhFxdCw4BHxA5r/KHlRC/I1M2uToNra+2vWActrnh8I1Lf9F9mmI/r61j4zs3Zp8U/164FDJT0JuAd4DfDaum3WAGfk/Q/HAI+U7V/o2Axuks4g6zxudndybw7iNTNLFEC1hTVDRIzn36XfBirA6oi4VdLp+fpPApeR3SN2B7AZeENCVh2bwW0/sqFVNwCrgW/nveaTXteKgpiZ9ZJocaiWiLiM7Mu/9rVP1jwO4K2zzON7M281syLDVf+XpL8DXkJWg31c0peBz0bEryPillYUpOU6GQArYVTNPZvKz5D2NLaUTgMwMlI+BlTKrHQAW8fL/zONqvyIle1Ke3+HO3h92/OX0glxmVJiOQFE9PbIrnqtvmLotDzU9j8DT6dmgFBEPLlI+kLvVl6T/S5fxoHFwFclnVe2wGZm/SBKLD3oc2Q3yI0DLyALaPr5oomLBNF7m6SfAOcBPwT+ICL+EvhD4E9TSmxm1tMiu2IouvSgeRFxFaCIuDsizgFeWDRxkevxJcArIuLu2hcjoirppFJFNTPrE63uY+iwrZKGgF/lnd730CQuXSNF+hiaRuSLiNuKZmRm1i+CbBKZPvYOYD7ZlJ7/QNacdGrRxL6PwcysgX6+YIiI6/OHm0gY9tpfQwXMzDqkn/sYJF3ZYIK1QsFPYYCvGMZv+m7pNFcs++OkvF61V/nhqvsuKn/qY2va27U94Zq4OnePpLwqCWM0NZEwNLZSakqQrujU98X2xG+muVH+czuReN4rCUOSuymi7/sYltRPsCapcB+DrxjMzBqollh6UFXSEyefSDqIEr9VBvaKwcxsNvr7goH3AD+QNHkn9POB04omdsVgZlYnu/O5f2uGiLginxr0mWQ34f+PiCg8R7KbkszMGujHO58lPTX/exTwRLKw3fcATywzh7SvGMzMGujF0UYFvJOsyeiDDdYFBe9+HtiKoXLUS0qneemifZLyigcfLZ3m/s3lR2nsM1Y+H4AYWVo6zdDWtLweSZj6ceFQ+dExqVJaB1KD4VXKTQ0JQCRcxJefJDbPKzEQYQolnItu68eWpIiY7Ec4ISK21q6TVHgOYjclmZnVCYJqiaUH/ajgaw0N7BWDmVmygIkeHYc6HUn7AQcA8yQdyc4L3j3IQmQU4orBzKxO0J9NScAfAa8nmy/6g+ysGB4Fziq6E1cMZmYN9GgT0bQi4kJJnwdOiYiLU/fjPgYzswaysBjFll4SEVXgLbPZR9uuGCStBk4CNkTE4flr5wBvBu7PNzsrnwe15YY23lU6zfcfKtwEN8VxE/eVTlNdvH/pNBrfVjoNwPqEaUSXTTyYlFdlzyVJ6TplKOFX4PbE308jHfrFOZ44rnK4gz8LU2JodVO/3+AGXCnpXcCXgMcnX4yIQv/Y7WxKugD4ONmUcrU+HBEfaGO+Zmaz06edzzXemP99a81rARSa87ltFUNEXCPp4Hbt38ysXfr9iiEinjSb9N3oYzhD0s2SVkta3GwjSadJWitp7f0bC4f4MDNrgWAiii+9RtJ8Sf9L0qr8+aFlpmLudMVwPnAIsAJYT+PbtgGIiFURsTIiVi5d0tvt1mY2WCavGIouPehzwBjw7Pz5OuB9RRN3tGKIiPsiYiLvNf80cHQn8zczKyTvYyi69KBDIuI8YDtARGyhRHSXjt7HIGlZRKzPn74cuKVtec0pP8LoqKUL0zJ7qGmLWFMjQwnDNCLtEzgnZfhJNe03w+aE6eKGqo/PvFGdsZG9SqcBGBkuf96HU94r0qJuKuG3WvLoooT4RT3667jl+r2PARiTNI/8YyjpEKDwsMZ2Dle9BDgWWCJpHXA2cKykFWSFvYtZjrU1M2uXXuw7KOEc4ApguaSLgeeQ3RFdSDtHJZ3S4OXPtis/M7NWya4Yul2KdBHxHUk/YedEPW8vM1GPQ2KYmdULmOjjmkHSGuASYE1ElG6rdUgMM7M6QfERST3aF/FB4HnAzyV9RdIry8zH4CsGM7MGJnry+76YiPge8D1JFbJZ294MrCYLvz2jga0YYstjpdOcf91vk/J6t24onebGpeXvzXgKabOq3TO+deaN6hw+VH6kEMDEnPJpQuUvXNMv88uPxNmSMNIKYH7CCKhqQvlSO0kTipccl6mSOLKrWwZgVBL5qKQ/Af4bcBRwYdG0A1sxmJkl6/8+hi8Bx5CNTPoEcHV+/1ghrhjMzOoMwBXD54DXRkT5yeVxxWBm1lCf9zFcMZv0rhjMzOp08opB0t5k8yYcTHbj76sj4qEG290FPAZMAOMRsbJdZfJwVTOzehFUq8WXWToTuCoiDgWuyp8384KIWNHOSgEG+YpheLR0klOOKD+rGkDlnmWl0yyZX7588eDm0mkAto+U/+BWE0Z1AWysjJVO89RK+by2ju5bOg3AHgnxpiYiMVZSQiyi7R1sv0gpX+oIqEpi5KhuCTralHQyWfggyEYOXQ28e7Y7lXQAcBA13/MRcU2RtINbMZiZzULJpqQlktbWPF8VEasKpt13MrhoRKyX9IQm2wXwHUkBfGq6/Ut6P9kw1Z+TNT1NpnfFYGaWIrtiKFUxbJyueUfSfwD7NVj1nhJ5PCci7s0rjisl3T7NFcDLgN+PiKSJ4l0xmJnVC1rRd7BzdxHHNVsn6b7JKQkkLQM2NNnHvfnfDZIuJZvPplnFcCcwQolQ27VcMZiZ1elwH8Ma4FTg3PzvN+s3kLQAGIqIx/LHLwH+fpp9bgZulHQVNZVDRLytSIFcMZiZNdDBG9zOBb4s6U3Ab4BXAUjaH/hMRJwI7AtcqmzAwDDwhRnuVViTL0lcMZiZ1YkIxjo0Z2dEPAC8qMHr9wIn5o/vBI4osc/CcZEaGdiKQaPzSqf5xm33JeX19tHyAec2DpUf1llZemDpNAAjY+WH/Wnv8kNwARgvn0QT5c/Fg1sSMgL2q5RPNx7lP0uQFnAuJT5PaoC61IB4KUYrKbdMde/W46DvYyUdCvwz8HRgR7jtiHhykfQDWzGYmaWKPg+iRxYr6Wzgw8ALgDdQ4sYQ3/lsZtbARDUKLz1oXkRcBSgi7o6Ic8jmZSjEVwxmZnWCnv3CL2qrpCHgV5LOAO4Bmt04twtXDGZm9fq/KekdwHzgbcA/kDUnnVo0cdsqBkmrgZOADRFxeP5aoSiCZmbd1O+dzxFxPYCkiIg3lE3fziuGC4CPAxfVvDYZRfBcSWfmz2cdLKpVXva0tMBsww83utN9essq5efAnHjg3tJpAObvfXD5RFvSphFdvEf5vCZGyp+/eZHYPZYQRK+SGMttJCHdRAfjxqVklVq8oS6OMErR753Pkp4FfBZYCDxR0hHAWyLir4qkb1vncx7D48G6l09m57yjF5LF8zAz6zl93vn8EeCPgAcAIuIm4PlFE3e6j6FoFEEknQacBrB8+fIOFc/MbCA6n4mI32pqaPXC03z27HDViFgVESsjYuXSJUu6XRwz241EZDcAFl160G8lPRsISaOS3gXcVjRxp68YCkURNDPrtj6/Yjgd+ChwALAO+A7w1qKJO10xzBhF0Mys2/q98zkiNgJ/lpq+ncNVLyGbrm6JpHVkt2c3jCLYDtW5i0qnuWXDpqS8nrxoYek0ew6XP/XxlKNLpwFYvGWkdJqJxYck5VXZVn7cSlTKn4s9K5XSaQDGRsp/LtieFkwtZerMlPhAcxIbhKsJY4xG0k47qpaPURVD3b3NKnUa014g6ULg7RHxcP58MfDBiHhjkfRtO/MRcUqTVbtEETQz6yUD0Pn8jMlKASAiHpJ0ZNHEvvPZzKxOvzclAUOSFk/eQJzfXFz4+94Vg5lZA31eMXwQ+JGkr+bPXwX8Y9HErhjMzOpkITE6M1FPO0TERZLWsjOi6isi4udF07tiMDOrF33fxwAwQhbFJPLHhQ1sxaDtW0unWbZoaVpmj20snWTzwvKzsVXW3146DcDw3keVT/PQuqS85u35tNJphrY8UjrNljkJo4uAfcYeSMgr7QZLVQvfaFqbKimvFNtTZotLLF63RxiV1e9B9CS9HXgz8DWyD9W/S1oVEf9aJH1/vVtmZh0weedzH3sTcExEPA4g6f3AjwFXDGZmKfr9ioHsKqH2knWCEpejrhjMzOr1/3DVzwHXSbo0f/4ysjDchbhiMDOr0+83uEXEhyRdDTyX7ErhDRHx06LpXTGYmTXQzxUDQETcANyQknZgK4bqvD1Lp1m/YVtSXhPLjyidZr/to6XTjO+zonQagAVj5YPpbF/6lKS8RsfLD1upJowwGq2kBQgKyp/3ocSROENjj5dOs4UFpdPMGU47F6MJcY/GJtK+LEcmyv9vxXD5WQ5bZQDufJ6Vga0YzMxmI1wxmJnZpIhgYqJ/73yeLVcMZmYN+IrBzMx2Cqi6YjAzs0kBxO7bkuSKwcyskejjGdxma2Arho2j5QPi/fF+DyXlNT66uHSaJ4yU/9BtjbS3a9Fo+bzGKvPS8koZAql9SqdZmPg/O6byw5jnJ45XHauUH4Y7P+G4to2n/bSdx/bSaaJckM4dNDFWPq8uDld1U5KZmdUJdz6bmdlOWR+DKwYzM5sUUHUfg5mZ1dqdrxjSgqzMkqS7JP1M0o35vKRmZj0lqlF4mQ1Jr5J0q6SqpJXTbHe8pF9IukPSmbPKdAbdvGJ4QUSUnxOzoK0JIzWqC8qPjgEYS7h1fm7CHImJ8csYTchrW2JmC4fKT2e5pVp+KNPcobTyPT5R/lykjkpKmQEsZaTQg+Np/8bDo+UDCpL4JdjVEUYJIqKTo5JuAV4BfKrZBpIqwCeAFwPrgOslrYmIn7ejQG5KMjNroFP3MUTEbQDStD9AjgbuiIg7822/CJwMtKVi6EpTElmn/3ck/UTSaY02kHSapLWS1t6/sW0XFmZmDUW1+AIsmfy+ypeG32uzcADw25rn6/LX2qJbVwzPiYh7JT0BuFLS7RFxTe0GEbEKWAXwh0cdtfv2AplZx0X5G9w2RsR0/QP/AezXYNV7IuKbBfbf6HKibd+LXakYIuLe/O+GfE7So4Frpk9lZtY5rRyVFBHHzXIX64DlNc8PBO6d5T6b6nhTkqQFkhZNPgZeQtb5YmbWG6Jzo5IKuh44VNKTJI0CrwHWtCuzblwx7Atcmne0DANfiIgrWp3J4rnlR7psTRyJsz0hVM38ifLTPjJcftpHgJFq+Tg12xJj4mxN+kiVP+/bIu03TUUJcaMSPxcpyTYmjDDaY07auRim/Ad3vGGLxsy0fUvpNFFJGDXVMtGxG9wkvRz4V2Ap8H8k3RgRfyRpf+AzEXFiRIxLOgP4NlABVkfEre0qU8crhrxXvfwkyWZmHdLJkBgRcSlwaYPX7wVOrHl+GXBZJ8rk4apmZvVi977z2RWDmVkDDrttZmZTeKIeMzPbIcLzMQykGW4vbyhxcAcJWRGUjx0zd6L8yA6AseHys7EtSBixkkn4Z0o5f0Npb1bKZGdbx9O+IBaMlC/jguHOfRk9Pl7+xKeNSYIYnpuYsnvclGRmZlNEtXxAyEHhisHMrF6EKwYzM9spcMVgZma1AmLCFYOZmU2KKtXx8qFkBsXAVgwPby1f2ydMdAbAPvPLn8ZqQiyioUgbKTSaECsJdS6+YgyVP3+pQ8xHEt7jOdVNaZmNlT+HUSn/uRgbSosptKi6uXxeI2nxurZe8oHSaea89j1JebWKm5LMzGwH9zGYmdlU4SsGMzObIqi6YjAzsx18H4OZmdXK5mNwxTBwUmbdOnCPtFnLKlsfLZ1GY+VHhJD4QY3R8rGStH1bUl4kjpwqa2juorSE1fHyeW15JCmroW3lZ+kb36PRfPHTm5twTABa/8vyee1zQFJeZ/z3fy+d5iPdHJUU4fsYzMxsKl8xmJnZTu5jMDOzqVwxmJlZjazzuTP9Zb3IFYOZWb3dvCmpcwFxakg6XtIvJN0h6cxulMHMbDpRnSi8DJqOXzFIqgCfAF4MrAOul7QmIn7eynyWL0o4tMTIbNU5CUMnU9J0UvkRrgOrOm9xt4vQHov27VhWH9l8W8fyaonwnc+ddjRwR0TcCSDpi8DJQEsrBjOzVIHnY+i0A4Df1jxfBxxTv5Gk04DTAJYvX96ZkpmZgfsYupBno4j4u7ThRMSqiFgZESuXLlnSgWKZmU0K9zF02Dqg9hLgQODeLpTDzKypQfzCL0qROhVWaobSMPBL4EXAPcD1wGsj4tZp0twP3N1k9RJgY6vLWVIvlAF6oxy9UAbojXL0QhmgN8rRyTIcFBFLZ7MDSVeQlbmojRFx/Gzy7CUdrxgAJJ0IfASoAKsj4h9nsa+1EbGyVWXr1zL0Sjl6oQy9Uo5eKEOvlKMXymDFdeUGt4i4DLisG3mbmdn0unKDm5mZ9a5BqBhWdbsA9EYZoDfK0QtlgN4oRy+UAXqjHL1QBiuoK30MZmbWuwbhisHMzFrIFYOZmU3RFxXDTNFYlflYvv5mSUe1oQzLJX1X0m2SbpX09gbbHCvpEUk35st721COuyT9LN//2gbrO3Eufr/mGG+U9Kikd9Rt05ZzIWm1pA2Sbql5bW9JV0r6Vf63YdS7VkX1bVKGf5F0e37OL5W0V5O0075/syzDOZLuqTnnJzZJ27Loxk3K8aWaMtwl6cYmaVtyLqwNIqKnF7J7HX4NPBkYBW4Cnl63zYnA5WThNp4JXNeGciwDjsofLyK7Sa++HMcC32rz+bgLWDLN+rafiwbvz+/Ibipq+7kAng8cBdxS89p5wJn54zOB96d8jmZZhpcAw/nj9zcqQ5H3b5ZlOAd4V4H3qyXnoVk56tZ/EHhvO8+Fl9Yv/XDFsCMaa0SMAZPRWGudDFwUmWuBvSQta2UhImJ9RNyQP34MuI0sIGCvafu5qPMi4NcR0ezO9JaKiGuAB+tePhm4MH98IfCyBkmLfI6SyxAR34mI8fzptWShXtqmyXkoomXnYaZySBLwauCS1P1bd/RDxdAoGmv9F3KRbVpG0sHAkcB1DVY/S9JNki6XdFgbsg/gO5J+kkegrdfRcwG8hub/+O0+F5P2jYj1kFXgwBMabNPJ8/JGsqu2RmZ6/2brjLw5a3WTJrVOnofnAfdFxK+arG/3ubBE/VAxFInGWihiaytIWgh8DXhHRDxat/oGsiaVI4B/Bb7RhiI8JyKOAk4A3irp+fVFbJCmXediFHgp8JUGqztxLsroyHmR9B5gHLi4ySYzvX+zcT5wCLACWE/WjLNLERu81q4x66cw/dVCO8+FzUI/VAxForF2JGKrpBGySuHiiPh6/fqIeDQiNuWPLwNGJLU0ZnhE3Jv/3QBcStY0UKuT0WtPAG6IiPsalLPt56LGfZPNZfnfDQ22aft5kXQqcBLwZxHR8Mu2wPuXLCLui4iJiKgCn26y7079rwwDrwC+1Gybdp4Lm51+qBiuBw6V9KT8F+prgDV126wB/iIfkfNM4JHJpoVWydtLPwvcFhEfarLNfvl2SDqa7Pw+0MIyLJC0aPIxWYfnLXWbtf1c1Gj6i7Dd56LOGuDU/PGpwDcbbFPkc5RM0vHAu4GXRsTmJtsUef9mU4bavqSXN9l3W89DjeOA2yNiXaOV7T4XNkvd7v0uspCNtPkl2WiK9+SvnQ6cnj8W2TzSvwZ+BqxsQxmeS3bJfTNwY76cWFeOM4BbyUZ6XAs8u8VleHK+75vyfLpyLvJ85pN90e9Z81rbzwVZRbQe2E726/dNwD7AVcCv8r9759vuD1w23eeohWW4g6ztfvKz8cn6MjR7/1pYhs/n7/nNZF/2y9p5HpqVI3/9gsnPQs22bTkXXlq/OCSGmZlN0Q9NSWZm1kGuGMzMbApXDGZmNoUrBjMzm8IVg5mZTeGKwczMpnDFYGZmU7hisJ4l6b/kAeHm5nfK3irp8Cbb/m0e2/8mSed2uqxmg8Q3uFlPk/Q+YC4wD1gXEf/cYJsTgL8DjouIzZL2joiUkNRmhisG63F5PJ/rga1kYTUmGmzzQbK4PJ/udPnMBpGbkqzX7Q0sJJs1b26TbUT7Qkeb7XZcMVivW0XWTHQx2ZSZjXwHeKOk+ZDNAd2hspkNpOFuF8CsGUl/AYxHxBckVYAfSXphRPxn7XYRcYWkFcBaSWPAZcBZnS+x2WBwH4OZmU3hpiQzM5vCTUnWNyT9AdlkNLW2RcQx3SiP2aByU5KZmU3hpiQzM5vCFYOZmU3hisGsz0h6uaTfStok6chul8cGjyuGGpLuknRcG/a7TNIaSfdKCkkH162fI2m1pEcl/U7SO+vWh6TH8y+CTZI+U7f+f+TpHsn3M2eG8jTdviaPyWVC0r9Os6+9JV2al+9uSa+tWTcq6av5eQ1Jx85QrpnOwwpJP5G0Of+7YhbH2bTcTfb1Ikm353l/V9JBNesk6f2SHsiX8yRpuv3N0geAMyJiYUT8tFU7lbRK0i8kVSW9fprt/jN/P5sOXpnpM2s9LiK85AtwF1kgtlbvd1/gr4BnkYVuOLhu/T8D3wcWA08DfgccX7M+gKc02fcfAfcBh+XprwbOnaYshbcHFgCbgOdPs79LgC+Rha14LvAIcFi+bhR4R/76euDYGc5T0/OQ7+tu4H8Ac4C35c9HU45zunI32NeSfP2ryMJy/Atwbc36twC/AA4EDgB+Dpzexs/peLPPwyz3+1bgRcBa4PVNtvkz4Jr8Mzk8zb6afma99P7S9QL0ykI2DLIKbMm/DP+2DXkMN6kY7gFeUvP8H4Av1jyfrmL4AvBPNc9fBPxumjIU3h44FbiTfPRag/ULgDHg9+rO4y4VDbCuQMXQ9DwAL8nXq2b9b6ipQIseZ5ly5+tOA35Ud9xbgKfmz38EnFaz/k21FccMx3xT/nmbXKLZeSKrECe3eRz4dZv+F37QqGIA9gR+CTzTFcNgL25KykXE68i+aP4kskv08+q3kfRESQ9Ps0zbHNGIpMXA/mRfEJNuIvulW+uavFnk63VNUYc1SLuvpH2aZFlm+1OBiyL/T2/g94CJiPjlDGWfUYHzcBhwc11Zbp5cL+m5kh6uWTfdcc5Y7vz9fG6jfUXE48Cv68o20/vXUEQckX/eFgLvJLvyuKHJttvy7QCOiIhDGm2nbA6LZp/RfytSrib+CTif7EquiGaf2Y7JmxA3SLqlRfubkHRjvqxpxT57kW9wKyEifgPs1eLdTv6jP1Lz2iNk0UQn/VfgWmA+8D7gW5JWRMR4nr4+LXn6B5rkN+P2kp6Y5/umGcr+SN1r9WUvaqbzMG1eEfEDpr430x3njOWOiPp93T/N9o3yWihJ01SqU+SV0PuA50bEo0XSNBMRz5hN+kYkrQSeA7ydrMlsJtN9ZjvpAuDjwEUt2t+WiFjRon31LF8xdN+m/O8eNa/tATw2+SQiromIsYh4mOwf80lkbfCT6evTAjwm6c9qOv8un2n7unL9BfCDiPi/ky9Iurxmf3/WYF+7lL2Emc5D2bymO87Z7mumsu0BbCpRKSwHvgycWncV0xMkDQH/Bry96Bf7DJ/ZjomIa4ApkzZJOkTSFfkAhu9Lemqny9XrXDFMNe0/ct6UVD9qZ1Pdl2W5DCMeIuuYPaLm5SOAW2co5+Sol1sbpL0vIh6IiIsnmyki4oSZtq/L4y+AC+vKekLN/i4ma28elnRoibI3PqCZz8OtwDPqRvs8Y5q8pjvOsuWesi9JC4BD6spW5v3bQdI84BvARyLi8hk2L0TZFKjNPqOfTNjlHsBK4EuSfkc2cRLAOknPK7iP2s9st60C/joi/hB4F1mlV9RcSWslXSvpZW0pXS/odidHLy1kl76ntWnfc8k6LQP4fWBuzbpzge+RjZ55KtkX5ORonMOAFUCFrMniI2Tt0CP5+uPJ2nyfnqf/T6YflTTj9sCzyTo3FxU4ri+SjfBZQNbUMGV0D1mH6VyyzueX5I+bdWZPdx4mRyW9Pd/nGUw/Kmna45yp3HX7Wpqv/9O8/O9n6qik04HbyEYk7U9WKZxes/5q4Jwm+74EuLjkZ6ktHbv5OZ4L/BB4c/54iOwLfb+a5b/kZTig0fmf6TPb6QU4GLglf7yQbODAjTXLbfm6VwC3NFi+XbOv/fO/TyYbxXhIN46p7ees2wXopQU4mawD+mHgXS3ed9QvNevmAKuBR8mGWL6zZt0L83+qx4ENZL8uD63b9zvzdI8CnwPmzFCWabcHPgV8vuBx7Z2X6fH83L22bv1dDY794Cb7anoe8vVHAj/J/7FvAI6sWfc8suabQsdZoNybgOfVPD8OuD3P++raYyD74jyPrMniwfxx7eipXwMvnuZzsZmpI5Oe12jbujTtqBiubvBeHdtgu4OpG5UEXA6cVfQz28mFqRXDHsD6Fu33AuCV3Tqudi4OomfWRpIOBL4SEc/qdll2V/mIqG9FxOH58x8BH46Ir+RNk8+IiJum20eebjGwOSK2SVoC/Bg4OSJ+3sbid4UrBjMbWJIuAY4lu0nxPuBssmbF84FlwAjZvTJ/X2Bfzya7mq6SNbF9JCI+256Sd5crBrMekw9i+FSDVXdHROl7RMzKcsVgZmZT9MUNbkuWLImDnvjEbhfDzPrADT/96caIWDqbfSzXvNhKtfD2Gxn7dkQcP5s8e0lfVAwHPfGJ/PCHP+x2McysD8ybP//u2e5jG1VerWWFt/+3uHvJbPPsJX1RMZiZdZKASpnI6QPWIu+KwcysgUqv3KfdBa4YzMzqlL5iGDCuGMzM6slXDAOp8mjRkPE7VecvTsssio9e6BtKjK+Yci5S8ko9553MK0VC+VRNi2QdlZHyeW1NCZwLY/ObTQ/S3HAXv5h9xWBmZlNkFUO3S9E9rhjMzHYhXzGYmdlOYveerMYVg5lZA75iMDOzHeRRSWZmVs9XDAOourB86JLtSjsdKcPqqgnT3w4l3ncfu/EHvN+kBDtO/VxsnSifbt7o/KS8Nm8vP+R3j9HutfJ7VJKZmU3h+xjMzGwXvmIwM7Mdss7n3bdmcMVgZtaArxjMzGwHdz4PqupE6SQjqXHjEu6RrCQEPouhtLdrqDqAQf46KBIDCqYEt4vKaOk0Q2OPl04DMGfOwtJptDUtYN+80XlJ6bpFDolhZma1JBgdcsVgZmY5NyWZmdku3JRkZmY7+IrBzMx24SuGNpA0F7gGmJPn89WIOFvSOcCbgfvzTc+KiMvaUIDSSaoJUx0CKCHATcrok1Spo2p6WWr8p6T3qoNfEOMJYY9Ghuck5bVtvPxotfmJn6UtCbGSRroZK0kw5IqhLbYBL4yITZJGgB9Iujxf9+GI+EAb8zYzmwWh3bgtqW0VQ0QEsCl/OpIvaWEgzcw6STC0G1cMbb1Wk1SRdCOwAbgyIq7LV50h6WZJqyUtbpL2NElrJa29f+PGdhbTzGwKAaoMFV4GTVuPKCImImIFcCBwtKTDgfOBQ4AVwHrgg03SroqIlRGxcumS8nMrmJklE6iiwsug6UhVFxEPA1cDx0fEfXmFUQU+DRzdiTKYmRUmMVQpvgyado5KWgpsj4iHJc0DjgPeL2lZRKzPN3s5cEtbCpAw+kQJ8ZWsO5T4myYlfhGJMapS8qqkjDAa214+DTA8XP64IhLPRfTfl6eGWve7WdJq4CRgQ0Qc3mC9gI8CJwKbgddHxA0tK0BJ7RyVtAy4UFKF7MrkyxHxLUmfl7SCrCP6LuAtbSyDmVlpan3n8wXAx4GLmqw/ATg0X44ha3I/ppUFKKOdo5JuBo5s8Prr2pWnmVmrtLLvICKukXTwNJucDFyUj+a8VtJeda0rHeU7n83M6kllRxstkbS25vmqiFhVIv0BwG9rnq/LX3PFYGbWC0TppqSNEbFyllnW69p9X64YzMzqCdTZ+RjWActrnh8I3NvJAtRyxVAjhipJ6VLi7ySJxJnYHCtpp4RzkRyXKSWvlI9SYoyvJImfwUg6sG6OZBJDnb1xbQ3Zjb9fJOt0fqRb/QvgisHMbFdqbeezpEuAY8n6ItYBZ5OFCSIiPglcRjZU9Q6y4apvaFnmCVwxmJnVUYsrhog4ZYb1Aby1ZRnOkisGM7MGOtyU1FNcMZiZ1dNgxkAqyhWDmVkdAUOdHZXUU1wx1HCspBqdHAGVkFdyrKSEvIK00WopkgZAJb9XaclSJJawe8RAhtMuyhWDmVkDgxg1tShXDGZm9dzHYGZmteSmJDMzq+emJDMz20ESlZHODTroNa4YzMzqtfjO534zuBVDyhC+Xg82lzosMUWvn4tE0ePHlRZrrrePqS+5j8HMzKZSS+d87jeuGMzM6mRzPrtiMDOzHUpP7TlQXDGYmdXbzfsY2nbkkuZK+v8k3STpVkn/O399b0lXSvpV/ndxu8pgZpYm62Mougyadl4xbANeGBGbJI0AP5B0OfAK4KqIOFfSmcCZwLtbnXlasLREHgG1U8pxdfJcdDCvlBFQqTOW9rpOzX7bMgJVdt/7GNr2XxKZTfnTkXwJ4GTgwvz1C4GXtasMZmYplPcxFF0GTVuPSFJF0o3ABuDKiLgO2Hdykuv87xOapD1N0lpJa+/fuLGdxTQzm0owNDRUeBk0bT2iiJiIiBXAgcDRkg4vkXZVRKyMiJVLlyxpWxnNzBrZna8YOjIqKSIelnQ1cDxwn6RlEbFe0jKyqwkzs96h3Xu4ajtHJS2VtFf+eB5wHHA7sAY4Nd/sVOCb7SqDmVkKgUcltcky4EJJFbIK6MsR8S1JPwa+LOlNwG+AV7WxDJ3RqZEuAzp6p9cpdUhNh85havynoZQhUIl5VftuVNLufcXQtoohIm4Gjmzw+gPAi9qVr5nZrO3mN7j5zmczswYcK8nMzHaQHF3VzMzquCnJzMx2cufzYIqh8ocWiYFqUgatDKVHZho8gxhrirTPU/IIKGs5NyWZmdkOkhjajYPouWIwM2vATUlmZraT72MwM7OpPFzVzMxqyFcMgymG55ROkx4TJ2E0U7XHR+J0cga3BKkjyFSdKJ/XUFonZMdGGCWe86RTOKhTzNXzcFUzM5tCYmh4pNul6Jrdt0o0M5vOUKX4UoCk4yX9QtId+Xz39euPlfSIpBvz5b0tP6aCfMVgZrYLQQs7n/PpBz4BvBhYB1wvaU1E/Lxu0+9HxEktyziRKwYzs3oCtfYGt6OBOyLiTgBJXwROBuorhp7gpiQzs12obFPSEklra5bT6nZ4APDbmufr8tfqPUvSTZIul3RYmw5uRr5iMDOrJwr3HeQ2RsTKGfZYr37Y2g3AQRGxSdKJwDeAQ8sUolUGtmLYOlF+qODcobThhUH5S86UQX+dHKLpaURboMeDA3YyXt9Qn41yVetvcFsHLK95fiBwb+0GEfFozePLJP2bpCURsbGVBSliQP8jzcxmYfKKoXWjkq4HDpX0JEmjwGuANVOylPaTsl9/ko4m+35+oLUHVszAXjGYmaVT2aakaUXEuKQzgG8DFWB1RNwq6fR8/SeBVwJ/KWkc2AK8JqI7cdhdMZiZNdDqWEkRcRlwWd1rn6x5/HHg4y3NNJErBjOzemrtFUO/aVsfg6Tlkr4r6TZJt0p6e/76OZLuqbm778R2lcHMLE3p4aoDpZ1XDOPA30TEDZIWAT+RdGW+7sMR8YE25s388U2l01RHFyTlVU1oBqxMbC+fj9I+gEPV8fJ5JQQhBBga31Y6zfah0dJphiNhpBVpAfHSgyv2eEDBFpdjoLT+Bre+0raKISLWA+vzx49Juo3GN3SYmfWY1obE6DcdOXJJBwNHAtflL50h6WZJqyUtbpLmtMm7CO/f2PFhvGa2O2v9cNW+0vaKQdJC4GvAO/IbOM4HDgFWkF1RfLBRuohYFRErI2Ll0iVL2l1MM7MaQkOVwsugaeuoJEkjZJXCxRHxdYCIuK9m/aeBb7WzDGZmSXbjpqS2VQz5HXyfBW6LiA/VvL4s738AeDlwS7vKYGaWRBrIK4Gi2nnF8BzgdcDPJN2Yv3YWcIqkFWQBpO4C3tKOzCuPrJ95ozrblv5eWl4Jo0KiUn52qJTRTwCVhNEx49W0vMqPL9o1klgRSpx6tBopo5J6e5rTftB/M4Lu3vcxtHNU0g9oPCLusgavmZn1DuGmJDMz20nI9zGYmVmN8vMxDBRXDGZmu3Afg5mZ1Wl1dNV+MrgVQwff1OGx8nGZIiEW0VDqL5hq+VE1lV6fcitxpFDS6JjUiPgJZUyJ5dRRiSOtJhIHdnXNbh5ddXArBjOz2diNhxu7YjAz24VcMZiZ2VThisHMzHZwH4OZme3Co5IM4PHtaUMnRofKn8btKp+m0uMDhYCkdtmkw0q8zE8JN5V62js6W1yC1JnfUlT67Ds2cFOSmZnVkjufzcysnisGMzPbyVcMZmZWx30MZmY2lSuGwROV8nOJ7TWa9kHQtrHSaSaGysdKGt36aOk0AFtGFpVOM7e6PSmvFClxmSI699Ht+dniOvgFVk2YeRCA8YnWFqTdpH6cdq5lBrZiMDObFV8xmJlZLfcxmJlZDfnOZzMzqyHclGRmZrV27/sY2nbkkpZL+q6k2yTdKunt+et7S7pS0q/yv4vbVQYzs2QaKr4MmBmvGCT9E3BeRDycP18M/E1E/K8Zko7n290gaRHwE0lXAq8HroqIcyWdCZwJvHsWx9C43A+uK53me4+l1VHHzV1fOs0DC+aVTrN84vHSaQBieGHpNJrYlpTX2MiC0mmGOxg4rpMjECcSjquSMFy1mhgeOiVgn6rjSXkNte83aNvszp3PRY78hMlKASAiHgJOnClRRKyPiBvyx48BtwEHACcDF+abXQi8rFyRzczabDKI3m56xVDkiCqSdtyNJWkeUOruLEkHA0cC1wH7RsR6yCoP4AlN0pwmaa2ktfdv3FgmOzOz2Zu8ya3IUmh3Ol7SLyTdkbeW1K+XpI/l62+WdFTLj6mgIhXDvwNXSXqTpDcCV7LzF/+MJC0Evga8IyIK37obEasiYmVErFy6ZEnRZGZmLdDaKwZJFeATwAnA04FTJD29brMTgEPz5TTg/NYeU3Ez9jFExHmSbgaOIxvE9Q8R8e0iO5c0QlYpXBwRX89fvk/SsohYL2kZsCGx7GZmbdPiPoajgTsi4k4ASV8ka1b/ec02JwMXRUQA10raa/K7spUFKaLQcNWIuAK4otE6ST+OiGc1eF3AZ4HbIuJDNavWAKcC5+Z/v1m20GZmbVeuYlgiaW3N81URsarm+QHAb2uerwOOqdtHo20OAHqzYpjB3CavPwd4HfAzSTfmr51FViF8WdKbgN8Ar2pBGXYRi/cvneaYPcuP3gGIR5udgub2mlN+JElsSwtg9vj28qNP5kf5wIAAkTAqKWWkSyRMp5qq10enpA7qShmgNTGU9hms0l9B9EIqO/XpxohYOc36Rjurf+eKbNMRrfjvaljwiPgBzT97L2pBvmZm7RHpFW4T64DlNc8PBO5N2KYjevunkJlZVwTVKL4UcD1wqKQnSRoFXkPWrF5rDfAX+eikZwKPdKN/AYrd4HYGWefxQ802aW2RzMy6r5UXDBExnn+XfhuoAKsj4lZJp+frPwlcRnaP2B3AZuANLSxCKUWakvYDrpd0A7Aa+Hbeaz7pdW0pmZlZlwRQbXHrfkRcRvblX/vaJ2seB/DW1uaaZsampDz0xaFkI4xeD/xK0j9JOiRff0tbS2hm1gURUXgZNEWHq4ak3wG/I4uBtBj4qqQrI+Jv21nAZB0ctULCqJp7NpWfOvNpbCmdBmBkpHwMqJTpSgG2jpf/JxlV+REr25X2/g53sOGz59tYE+IypcRyAojor+7Mdlwx9JMifQxvI7vfYCPwGeB/RsR2SUPAr4DerBjMzFIFTLhimNYS4BURcXftixFRlXRSe4plZtZdg9hEVFSRkBjvnWbdba0tjplZ9wWQ1mg2GDyDm5lZA7vxBYMrBjOzRtz5PIDGb/pu6TRXLPvjpLxetVf5UUn7Lip/6mNr2tu1PeGauDp3j6S8KglDcTSRMAKqUmpKkK7o1PfK9sRvsLlR/nM7kXjeKwkjz7opwn0MZmZWx30MZmY2xW58weCKwcysXnaD2+5bM7hiMDNrYPetFlwxmJk15FFJA6hy1EtKp3npon2S8ooHHy2d5v7N5Udp7DNWPh+AGFlaOs3Q1rS8HkmY4WvhUPnRMalSWgdSYx5Vys0ABkAkTJFSfi7APK/EeFMplHAuum03bkka3IrBzCxVEFR348YkVwxmZvVaP7VnX3HFYGbWgPsYzMxsh8BXDGZmVsd9DG0gaTVwErAhIg7PXzsHeDNwf77ZWfk8qC03tPGu0mm+/9D8pLyOm7ivdJrq4v1Lp9H4ttJpANYnzBa3bOLBpLwqey5JStcpQwn/7NsTRgoBjHToi2U8sc1juIOTqqXE0Oq23fmKoZ0fjQuA4xu8/uGIWJEvbakUzMxmY/LO56LLoGnbFUNEXCPp4Hbt38ysbQImduMoet2YofsMSTdLWi2p6Sz1kk6TtFbS2vs3buxk+cxsN7e7XzF0umI4HzgEWAGsBz7YbMOIWBURKyNi5dIlvd1ubWaDJpiI4sug6eiopIjY0Usr6dPAtzqZv5lZEY6u2kGSlkXE+vzpy4FbOpm/mVkhu3kfQzuHq14CHAsskbQOOBs4VtIKsgr5LuAtbct/Tvmhp0ctXZiW2UNNu0qaGhlKGL8XaZ/UOSnjEqtprYybE+YRHao+XjrN2MhepdMAjAyXP+/DKe8VaWGbldC6mzzsNCGw3e7yK9pXDG0SEac0ePmz7crPzKyVBrHvoCjf+WxmVie7Yuh2KbrHFYOZWZ2IYPtu3MngisHMrIEJXzGYmdkkdz4PqNjyWOk051/326S83q0bSqe5cWn5m/aeQtp0m/eMby2d5vCh8iOFACbmlE8TKj+sZiK5Abj8SJwtCSOtAOYnjICqJpQvtZM0oXjJAfsqiSO7uiZm8xnrfwNbMZiZpfIVg5mZ7cJ9DGZmtoOvGMzMbKoIqh3qY5C0N/Al4GCyiBCvjoiHGmx3F/AYMAGMR8TKdpWpG2G3zcx6WpA1JRVdZulM4KqIOBS4Kn/ezAvySc7aVinAIF8xDI+WTnLKEeWn2wSo3LOsdJol88uXLx7cXDoNwPaR8p/casKoLoCNlbHSaZ5aKZ/X1tF9S6cB2CMh3tREJMZKSohFtL2DDdsp5UsdAVVJjBzVTR1sSjqZLK4cwIXA1cC7O5V5I4NbMZiZJcquGEpVDEskra15vioiVhVMu+9k1OmIWC/pCdMU6zuSAvhUif2X5orBzKxeULaPYeN0zTuS/gPYr8Gq95TI4zkRcW9ecVwp6faIuKZMIYtyxWBmVmeyj6Fl+4s4rtk6SfdNzlUjaRmwock+7s3/bpB0KXA00JaKwZ3PZmYNdHDO5zXAqfnjU4Fv1m8gaYGkRZOPgZfQxonOfMVgZlYn6OhczucCX5b0JuA3wKsAJO0PfCYiTgT2BS5VNmBgGPhCRFzRrgINbMWg0Xml03zjtvtm3qiBt4+Wjyu0caj86J3K0gNLpwEYGSs/ukN7lx9pBcB4+SSaKH8uHtySkBGwX6V8uvEo/1mCtLhCKfF5UuMQpcY9SjFaSWmc6OINZh2MlRQRDwAvavD6vcCJ+eM7gSM6UiAGuGIwM0sVOIiemZnVCEdXNTOzeq4YzMxshyBcMZiZWQ03JbWHpNXAScCGiDg8f61QFMFuednT0uLvDD/c6IbG6S2rlJ/qbOKBe0unAZi/98HlE21Jmy1u8R7l85oYKX/+5kXiLTgJsZIqiSF7RhLSTXQwPFBKVqnFG+rmCKMEu3vncztvcLsAOL7utTJRBM3MumKy87noMmjaVjHkMTwerHv5ZLLogeR/X9au/M3MZmN3rhg63cdQNIogkk4DTgNYvnx5h4pnZubO556NlRQRqyJiZUSsXLpkSbeLY2a7kYjszvCiy6Dp9BVDoSiCZmbd5iuGzpkxiqCZWbft7p3P7RyuegnZdHVLJK0DzqZJFMF2qM5dVDrNLRs2JeX15EULS6fZc7j8qY+nHF06DcDiLSOl00wsPiQpr8q28gMao1L+XOxZqZROAzA2Uv5zwfbyQ1whberMlMBxcxJ/3lUTBp+OpJ12VC0fvDCGunebVRCMTaS974OgbWc+Ik5psmqXKIJmZr3EsZLMzGwXrhjMzGyH7M5nNyWZmdmkGMxO5aJcMZiZ1dndYyUNbMWg7VtLp1m2aGlaZo9tLJ1k88Ly03RW1t9eOg3A8N5HlU/z0LqkvObt+bTSaYa2PFI6zZY5CaOLgH3GHkjIK+0GS1UnUlIl5ZVie8o0oonF6+YIoxSTN7jtrvrr3TIz6wBfMZiZ2VQermpmZrV29yB6rhjMzBpwxWBmZjv4zucBVZ23Z+k06zdsS8prYvkRpdPst320dJrxfVaUTgOwYKx8MJ3tS5+SlNfoePlhK9WEEUajlbQAQUH58z6UOBJnaOzx0mm2sKB0mjnDaediNCHu0dhE2pflyET5/60YLj/9bSuFKwYzM5sUAVVXDGZmtlMQ4YrBzMxquCnJzMx2clOSmZnVCiB23+Cqg1sxbBwtH/foj/d7KCmv8dHFpdM8YaT8r5GtkfZ2LRotn9dYZV5aXikjXbRP6TQLE3/Mjan8aLX5icOSxirlR1vNTziubeNp32Dz2F46TVB+NkAATYyVz6vbo5Lcx2BmZju4KcnMzKYKdz6bmdlOWR+DKwYzM5sUUN2N+xjS7qU3MxtwUY3Cy2xIepWkWyVVJa2cZrvjJf1C0h2SzpxVpjPoyhWDpLuAx4AJYDwimp6MVFsTRmpUF5QfHQMwNlE+r7kJU2ElhqlhNCGvbYmZLRwqP2vZlmr5oUxzh9LK9/hE+XOROiopZQawlJFCD46n/RsPj5aPG0Xil2C3Rxil6GBT0i3AK4BPNdtAUgX4BPBiYB1wvaQ1EfHzdhSom01JL4iI8nNimpm1WUR0bFRSRNwGIE37A+Ro4I6IuDPf9ovAycDAVQxmZj2r5H0MSyStrXm+KiJWtbA4BwC/rXm+DjimhfufolsVQwDfkRTApxqdQEmnAacBLF++vMPFM7PdXck7nzdO1yQu6T+A/Rqsek9EfLPA/htdTrTtkqZbFcNzIuJeSU8ArpR0e0RcU7tBXlmsAvjDo47afYcHmFnHRcBE4h3ljfcXx81yF+uA2l/IBwL3znKfTXVlVFJE3Jv/3QBcStZ+ZmbWG6Jzo5IKuh44VNKTJI0CrwHWtCuzjl8xSFoADEXEY/njlwB/3+p8Fs8tP9Jla+JInO0JPyzmT5Sf3Yvh8rN7AYxUy8ep2ZYYE2dr0keq/HnfFmm/aSpKiBuV+LlISbYxYYTRHnPSzsUw5T+44w1bNGam7VtKp4lKwqiplomO3ccg6eXAvwJLgf8j6caI+CNJ+wOfiYgTI2Jc0hnAt4EKsDoibm1XmbrRlLQvcGneAz8MfCEiruhCOczMGurknc8RcSlZy0n96/cCJ9Y8vwy4rBNl6njFkA+3Kj9JsplZp4RDYpiZWR1HVzUzsyk8H4OZme0Q4bDbZmZWx01JA2iGuCMNJY76IyErgvJBxeZOlB/yBzA2XH6azgUJQxkzCf9MKedvKO3NSrlnaet42hfEgpHyZVww3Lkvo8fHy5/4tMGqEMNzE1N2T1TLB4QcFANbMZiZJYtwxWBmZjsFrhjMzKxWQEy4YjAzsx18xWBmZrXcxzCYHt5a/k1NmAETgH3mlz+N1YQgdUMlA8RPGk0Iooc6F3g3hsqfv9R7j0YS3uM51U1pmY2VP4dRKf+5GBtKCza3qLq5fF4jaYEct17ygdJp5rz2PUl5tYorBjMz28Gdz2ZmNlX4isHMzKYIqq4YzMxsB3c+m5lZrWyiHlcMAydlOsYD90ibzrKy9dHSaTRWfkQIiR/UGC0fK0nbtyXlReLIqbKG5i5KS1gdL5/XlkeSshraVn761vE99iudZm7CMQFo/S/L57XPAUl5nfHf/710mo90c1RShG9wMzOzqXzFYGZmO7mPwczMpnLFYGZmNbLO5870l/UiVwxmZvXclNR5ko4HPgpUgM9ExLmtzmP5ooRDSwzAU52TMEImJU0nlR/INLCq8xZ3uwjtsWjfjmX1kc23dSyvVnHF0EGSKsAngBcD64DrJa2JiJ93uixmZo1EBBPbE4JPDohuXDEcDdwREXcCSPoicDLgisHMeoObkjruAOC3Nc/XAcfUbyTpNOA0gOXLl3emZGZmwO4+KqlzQfd3ahQRf5fG/YhYFRErI2Ll0iVLOlAsM7OdojpReBk03bhiWAfUXgIcCNzbhXKYmTW2mzclKVKnwkrNUBoGfgm8CLgHuB54bUTcOk2a+4G7m6xeAmxsdTlL6oUyQG+UoxfKAL1Rjl4oA/RGOTpZhoMiYulsdiDpCrIyF7UxIo6fTZ69pOMVA4CkE4GPkA1XXR0R/ziLfa2NiJWtKlu/lqFXytELZeiVcvRCGXqlHL1QBiuuK/cxRMRlwGXdyNvMzKbXjc5nMzPrYYNQMazqdgHojTJAb5SjF8oAvVGOXigD9EY5eqEMVlBX+hjMzKx3DcIVg5mZtZArBjMzm6IvKgZJx0v6haQ7JJ3ZYL0kfSxff7Oko9pQhuWSvivpNkm3Snp7g22OlfSIpBvz5b1tKMddkn6W739tg/WdOBe/X3OMN0p6VNI76rZpy7mQtFrSBkm31Ly2t6QrJf0q/9swHOpMn6NZluFfJN2en/NLJe3VJO20798sy3COpHtqzvmJTdK25DxMU44v1ZThLkk3NknbknNhbRARPb2Q3evwa+DJwChwE/D0um1OBC4nC7fxTOC6NpRjGXBU/ngR2U169eU4FvhWm8/HXcCSada3/Vw0eH9+R3ZTUdvPBfB84CjglprXzgPOzB+fCbw/5XM0yzK8BBjOH7+/URmKvH+zLMM5wLsKvF8tOQ/NylG3/oPAe9t5Lry0fumHK4Yd0VgjYgyYjMZa62TgoshcC+wlaVkrCxER6yPihvzxY8BtZAEBe03bz0WdFwG/johmd6a3VERcAzxY9/LJwIX54wuBlzVIWuRzlFyGiPhORIznT68lC/XSNk3OQxEtOw8zlUOSgFcDl6Tu37qjHyqGRtFY67+Qi2zTMpIOBo4Ermuw+lmSbpJ0uaTD2pB9AN+R9JM8Am29jp4L4DU0/8dv97mYtG9ErIesAgee0GCbTp6XN5JdtTUy0/s3W2fkzVmrmzSpdfI8PA+4LyJ+1WR9u8+FJeqHiqFINNZCEVtbQdJC4GvAOyLi0brVN5A1qRwB/CvwjTYU4TkRcRRwAvBWSc+vL2KDNO06F6PAS4GvNFjdiXNRRkfOi6T3AOPAxU02men9m43zgUOAFcB6smacXYrY4LV2jVk/hemvFtp5LmwW+qFiKBKNtSMRWyWNkFUKF0fE1+vXR8SjEbEpf3wZMCKppTHDI+Le/O8G4FKypoFanYxeewJwQ0Tc16CcbT8XNe6bbC7L/25osE3bz4ukU4GTgD+LiIZftgXev2QRcV9ETEREFfh0k3136n9lGHgF8KVm27TzXNjs9EPFcD1wqKQn5b9QXwOsqdtmDfAX+YicZwKPTDYttEreXvpZ4LaI+FCTbfbLt0PS0WTn94EWlmGBpEWTj8k6PG+p26zt56JG01+E7T4XddYAp+aPTwW+2WCbIp+jZMrmMX838NKI2NxkmyLv32zKUNuX9PIm+27reahxHHB7RKxrtLLd58Jmqdu930UWspE2vyQbTfGe/LXTgdPzxyKbR/rXwM+AlW0ow3PJLrlvBm7MlxPrynEGcCvZSI9rgWe3uAxPzvd9U55PV85Fns98si/6PWtea/u5IKuI1gPbyX79vgnYB7gK+FX+d+982/2By6b7HLWwDHeQtd1PfjY+WV+GZu9fC8vw+fw9v5nsy35ZO89Ds3Lkr18w+Vmo2bYt58JL6xeHxDAzsyn6oSnJzMw6yBWDmZlN4YrBzMymcMVgZmZTuGIwM7MpXDGYmdkUrhjMzGwKVwzWsyT9lzwg3Nz8TtlbJR3eZNu/zWP73yTp3E6X1WyQ+AY362mS3gfMBeYB6yLinxtscwLwd8BxEbFZ0t4RkRKS2sxwxWA9Lo/ncz2wlSysxkSDbT5IFpfn050un9kgclOS9bq9gYVks+bNbbKNaF/oaLPdjisG63WryJqJLiabMrOR7wBvlDQfsjmgO1Q2s4E03O0CmDUj6S+A8Yj4gqQK8CNJL4yI/6zdLiKukLQCWCtpDLgMOKvzJTYbDO5jMDOzKdyUZGZmU7gpyfqGpD8gm4ym1raIOKYb5TEbVG5KMjOzKdyUZGZmU7hiMDOzKVwxmJnZFK4YzMxsiv8few+t8z9yeTAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(2, 1)\n", "ds.woce.isel({'z_f':1}).plot(ax=ax[0])\n", "w.isel({'z_f':1}).plot(ax=ax[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 2 fields look similar, which is confirmed by computing the difference." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray ()>\n",
       "array(2.07489191e-17)
" ], "text/plain": [ "\n", "array(2.07489191e-17)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs((w - ds.woce)).max()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Vorticity\n", "Here we compute more derived 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 v}{\\partial x} - \\frac{\\partial u}{\\partial y} \\ . $$\n", "\n", "The NEMO discretization is [1]\n", "\n", "$$\n", "\\zeta = \\frac{1}{e_{1f} e_{2f}} \\left[\\delta_{i+1/2}(v \\cdot e_{2v}) - \\delta_{j+1/2}(u \\cdot e_{1u}) \\right] \n", "$$\n", "\n", "
\n", "[1] NEMO book v4.0.1, pp 22\n", "
\n", "\n", "In xgcm, we calculate this quantity as" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Coordinates:\n", " * x_f (x_f) float64 0.5 1.5 2.5 3.5 4.5 5.5 ... 15.5 16.5 17.5 18.5 19.5\n", " * y_f (y_f) float64 0.5 1.5 2.5 3.5 4.5 5.5 ... 35.5 36.5 37.5 38.5 39.5\n", " * t (t) object 1050-07-01 00:00:00\n", " * z_c (z_c) float32 5.0 15.0 25.0 35.0 ... 3.38e+03 3.787e+03 4.22e+03" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "zeta = 1/(domcfg.e1f*domcfg.e2f) * (grid.diff(ds.vo*domcfg.e2v, 'X', **bd) - grid.diff(ds.uo*domcfg.e1u, 'Y', **bd)) * domcfg.fmaskutil\n", "zeta.coords" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\zeta$ is located in the F point (*x_f*, *y_f*).\n", "\n", "We plot the vertical integral of this quantity, i.e. the barotropic vorticity:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAI/CAYAAAClCKSgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4HUlEQVR4nO3df/RddX3n+9cr3wQkxAgkBL8lYFDj2OhYoLkMc72jVsQGbqex9WrhTsfg4l5GR7raTufegbFd17rqKqOdudU7KDcqLdgfyGCtcUpFTLVOO/4gCqKAQMQ0Br4SEsAYw4Xwzfv+cfaBw8n5sc+Pvfdn7/N8rJX1PWef/eNzvjlnv77vz2f/cEQIAIAiLKm6AQCA5iJkAACFIWQAAIUhZAAAhSFkAACFIWQAAIVZWubGVq1aHae96PQyNwmgpr51++37IuLkSdbx4iXL4wktTqtJA/0wnrolIjaVsrEaKTVkTnvR6fqbv/27MjcJoKZWrTz+HyZdxxNa1NuXrZ1Gc4b6/aceWF3KhmqG7jIAQGFyhYztE2zfZPu7tu+x/U9tn2T7Vtv3Zz9PLLqxAIB6yVvJfFDS5yLi5ZJ+RtI9kq6QtD0i1kvanj0HAOAZQ0PG9kpJr5H0cUmKiKci4nFJmyVdl812naQ3FdNEAEBd5alkXizpEUl/ZPt22x+zfbykUyJiQZKyn2sKbCcAoIbyhMxSSWdL+khEnCXpJxqha8z2ZbZ32N6xf9++MZsJAKijPIcw75G0JyK+lj2/Sa2Qedj2fEQs2J6XtLfXwhGxVdJWSTrz7LPHuq/AiiOHdHDJ8nEWHcuKI4dK2xYwrrK/E2VuD80xtJKJiB9K+oHtf5RNOk/S3ZK2SdqSTdsi6TPTbtyKI4ee2eGXteMnYFAXZX8nOr+PQF55T8b8NUl/avsYSQ9IertaAXWj7Usl7Zb0lmk1qqoPMl8g1E37M1tEldHv+0BVg1HkCpmIuEPSxh4vnTfV1mjwjr6oDzfhgrqYO7DwzOPFlfPPPJ72d6PzO9Frm0WGG5ql1MvKDJJ3R1/klwlITecOvt9rnTv+Sb8b3d+H7u0XsU00WxIhM+qOflofbAIG09C9453GukZdZhoVRr/qRZJi725JktecPtVtovkqDZlJdvKTBg0Bg3F26GWur5f2zr6tc6cvjV9h9AuY7u3F3t1T2yZmQ2UhM42d/LgfagImTWXspOuoe0ff/Vp7py+NXmEM6h7r3O7hhV1aNr/uOdOpapBH6SEz7R38JH+xYTh2/NUYFCxS/p2+1P87krd6Obyw6zk/O7c7qKqRCBuUHDJzcaSQ9eYNGgKG0EjdoHBp7+Q7n7d3+O1l83Zl5QmY7u312u64AYfZ0Zj7yQw7UYyAAYDyJXF02TT1+suJgHlWvyOgqHDS0K4IelU07eqhu9uqe9m2cY5285rTn9n2svl1fauZQdsFOjUuZKTnBg0Bk8+0w6fX+giy/Dp39t26w6U9f6fu3/+gLqvFlfMj/d/02j7QTyNDRiJcpmUa5370WxehM9igqqbXfG15AubgkuV9vyODqplh1RPQrbEhg/RNo9oZJQRTCLVxKsZBYTONLrL2cv3akLfbbJLto7kIGSSlyJ1UmeNRo76P9vyjhE2egBnUTZa3mpHyddFhMNub1LqV/Zykj0XEVV2vO3v9QkmHJF0SEd8ctKztkyR9UtI6SbskvTUiHrO9TNLH1LoX2FJJ10fE72fLXCzp30sKSQ9J+tWI2Gf7dLXucnxCtp0rIuLmSd83IYOZl9Jf36OETa/lOo166PCoYzPIz/acpKslna/WPbpus70tIu7umO0CSeuzf/9E0kck/ZMhy14haXtEXGX7iuz5v1PrqvjHRsQ/tr1c0t22/zxb/oOSNmTB8n5Jl0t6j6TflnRjRHzE9gZJN6sVXhNJ/hDmuQMLfPAxcxZXzucOv1FCcv8Ti7nnHVSpUMWM7BxJOyPigYh4StINkjZ3zbNZrYojIuKrkk7Ibgg5aNnNalUfyn6+KXscko63vVTScZKeknRAkrN/x2eV00q1qpn2Miuzxy/omD6RpCuZznCZ5kUI826zl5T+6kXz5alseulVxfQKmO4us0mqGb4bA50q6Qcdz/eoVa0Mm+fUIcueEhELkpTdpXhNNv0mtQJoQdJySb8ZEY9Kku13Svq2pJ9Iul/Su7Jl3iPp87Z/TdLxkt4wzhvtlmwlk2r1kmq70Gz9KptpdJMN0qtioYrpa7XtHR3/Lut4zT3m774dfb958izb7RxJi5J+StIZkn7L9ouzsZp3Sjore+1OSVdmy1ws6Y8jYq1a40KfsD1xRiRdyfTSfdmKMvGXGqrWWWn0+zz2u5TLquPmRuou66VuAfO8JUv0shXHlLOxR7UvInrd3FFqVR+ndTxfq6O7o/rNc8yAZR+2PZ9VMfOS9mbT/1dJn4uIw5L22v57tW48uUqSIuJ7kmT7RrXGcSTpUkmbste/Yvt5klZ3rHMsyVYyqWj/BUnAIBWTfBZXHTc38vrzBAvfj6Fuk7Te9hnZbewvkrSta55tkt7mlnMl/SjrChu07DZJW7LHWyR9Jnu8W9Lrs3UdL+lcSd+V9KCkDbZPzuY7X9I9HcucJ0m2f1rS8yQ9MukbT7KSGdYlVWQ10/5LkS8NUlb257NuFUxqIuJp25dLukWtw4OvjYi7bL8je/0atY7mulDSTrUOYX77oGWzVV8l6Ubbl6oVEm/Jpl8t6Y8kfUet7rY/iog7Jcn270r6su3Dkv5B0iXZMr8l6aO2f1Ot7rhLImJYt9xQSYZM1QgYzIppXBmD70s+2TknN3dNu6bjcejZQfihy2bT9yurPrqmH9SzgdP92jWSrukx/W5Jrx74JsZQ2+4yBuBRd1wCH7MguZAhPDAL2gGTatAUcY4OZlNyIQM02cEly48KllSDBpiGpEJm1CqGqgd10StcipRnrGXS8RiqGOSRVMgATZQnXKhm0FSEDBqnXTWUXT1Mooj7HxVZqVDFIC9CBo2Wws3rVhw5VFk7+m03b3vyXsoG6IeQQeO0d6ApBEynSXf4097uOAgYjIqQQSOlFjBtVVU1ndscdfvtYCFgMA5CBqhAFdXWJAFHwGBchAxQkVSrLWCaCBkAQGEIGQBAYQgZAEBhah8yXFoGANJV+5ABAKSLkAEAFIaQAQAUJpmQYWwFAJonmZABADQPIQMAKAwhAwAoDCEDACjM0qobAABFWbZsiU5dU9LdUR8tZzN1QyUDACgMIQMAKAwhAwAoDCEDACgMIQMAKAwhAwAoDCEDACgMIQMAKAwhAwAoDCEDACgMIQMAKAwhAwAoDCEDACgMIQMAKAwhAwAoDCEDACgMIQMAKAwhAwAoDCEDACgMIQMAKAwhAwAoDCEDACgMIQMAKAwhAwAoDCEDACgMIQMAKMzSqhsAAEVZeuycTnrpSeVs7LvlbKZuqGQAAIUhZAAAhSFkAACFIWQAAIUhZAAAhSFkAACFIWQAAIUhZAAAhcl1MqbtXZJ+LGlR0tMRsdH2SZI+KWmdpF2S3hoRjxXTTABAHY1SyfxcRJwZERuz51dI2h4R6yVtz54DAPCMSbrLNku6Lnt8naQ3TdwaAECj5A2ZkPR529+wfVk27ZSIWJCk7OeaIhoIAKivvBfIfHVEPGR7jaRbbee+FFwWSpdJ0umnre05z9yBhYHriL275TWn93197sCCFlfO520SAKAkuSqZiHgo+7lX0qclnSPpYdvzkpT93Ntn2a0RsTEiNq5etWrkBsbe3SMvAwBIw9CQsX287ee3H0t6o6TvSNomaUs22xZJnymqkQCAesrTXXaKpE/bbs//ZxHxOdu3SbrR9qWSdkt6S3HNBADU0dCQiYgHJP1Mj+n7JZ1XRKMAAM3AGf8AgMIQMgCAwtQmZDjKDADqpzYhAwCon6RDZpTqZdgJnQCA8iUdMgDQFLY32b7X9k7bR11Q2C0fyl6/0/bZw5a1fZLtW23fn/08MZu+zvYTtu/I/l3TscyvZOu/y/b7O6b/G9t3Z69tt/2iabxvQgYACmZ7TtLVki6QtEHSxbY3dM12gaT12b/LJH0kx7KDrob/vezK+WdGxDuyda2S9AFJ50XEKySdYrt9KsrtkjZGxKsk3STp/ZqCWoUMg/8AauocSTsj4oGIeErSDWpdyb7TZknXR8tXJZ2QXbJr0LKjXg3/xZLui4hHsudfkPRmSYqIL0bEoWz6VyX1vtjkiGoVMgBQU6dK+kHH8z3ZtDzzDFp20NXwz7B9u+2/tf3Psmk7Jb08605bqlYondajvZdK+uuc722gvFdhLky/AXuqFgCTmjt2mU586cllbW617R0dz7dGxNbssXvMH13P+82TZ9luC5JOj4j9tn9W0l/afkVEPGb7nWrd1fiIpP+uVnXzbCPsX5W0UdJrh2wjl8pDZpq45D+ACu3ruHNwtz16bsWwVtJDOec5ZsCyD9uej4iFzqvhR8STkp7MHn/D9vckvUzSjoj4rKTPSs/cimWxvWLbb5D0bkmvzdYxsdp1l1HhAKih2yStt32G7WMkXaTWlew7bZP0tuwos3Ml/SjrAhu0bM+r4ds+OTtgQLZfrNbBBA9kz9dkP0+U9K8lfSx7fpak/1fSL2a3dZmKRlUys+jgkuVaceTQ8BkBVCYinrZ9uaRbJM1JujYi7rL9juz1ayTdLOlCtcZNDkl6+6Bls1Vfpd5Xw3+NpPfaflqtSuUdEfFo9toHbbcvevzeiLgve/wBSSsk/Zfsqvu7I+IXJ33vhEyNHVyy/JmfBA2Qtoi4Wa0g6Zx2TcfjkPSuvMtm03teDT8iPiXpU33WdXGf6W8Y0Pyx1a67DC3tgOn3HABSQMjUUL9AObhkOWEDICmETM3kCRGCBkAqkgwZjiDrbZTwIGgApCDJkMHRxgkNus8AVK2WITNrlc6kQUHQAKhKLUNmlkwrIAgaAFUgZBI27WCg+wz8/6NsnIyZKHYGmIZen6P2NE7gRRkImQQVHTBcIaC5Rj0Ckc8BikbIJKaMCoYdS/OM+7kpuqqhIgchkxC+kKjKNMOGzzE6ETKJ4IuJcU3zszNO2PDZxSCETAL4kiI1w8KGzyzyImRmDOMxzVLGQSJS63NDsGAclZ4nM3dgocrNJ4EvLsZV5meHzynGldzJmLN2yRgAaLLkQgbAcFQWqIvahkwTKh52FACarrYhg9Ex6N8M/HGCOmlcyNTlYAJ2FABmQeNCBmgy/jhB3XCeTAXYUWBc3V2efJaQOkIGqDHOyB9s7thlesFLTq26GTONkAEaqFf4EDyoAiFTMr7oqEpn8PA5RFkY+AdmEIezoyyEzIxgp4JufCZQBkKmRHRRIDUEDYpGyAAzjqBBkQgZAAQNCkPIlISuMqSOoEERCBkAzyBoMG2EDIDnIGgwTYRMCegqQ90QNJgWQgYAUJhKQ2Zx5fxR07zm9ApaAgAoQiMrmbrcuKwsdH1gHHxuMA2NDBkAQBoIGQB9Uc1gUrUNGcZuACB9tQ0ZAED6CBkAA9FlhkkQMgCAwhAyAIaimsG4CBkAQGGWVt2AXrzmdMXe3VU3oxH4CxTTMupniWv2QaKSAVCQFUcO8UcO0qxkmoS/5jAr+KyjF0KmQHzp0DR8pjGqWoZMHc7258uIJqnr53nJscdo+bp1VTdjpjEmU4C6fiGBXvg8YxK1rGRSxhcSTcFnGdNQeSXT68ZlUj26xLrxpURT8FnGtNSukqlj+AB1Qbhg2moXMqniyzmacX9fnHdRHD7DKAIhMwV8OVvK+D10b6Oq0Bn0XusYhHyGURRCBiNJbWdUdOiM835TCcI8Uvv/RPPUKmRSHI+ZlS9pXd5niu1MNXRS/F2heWoVMqgGO6Ppqjp0+P9EmZIOmdSvxjwLX9ZZeI9V6/wdFxE4/B+iSkmHTMpm4Ys7C+8xNe3f+aRhw/8dUlH5yZiDpFrF1OkLPG5b6/Qem+jgkuVjH3TA/x1SknTIpGgWvsCz8B7rYpT/C/7fkKJahUzVlU1dv8TsqOotT3XC/xtSlTtkbM/Zvt32f82en2T7Vtv3Zz9PLK6Z1ZuFL/EsvMc66/X/Q/cYUjdKJfPrku7peH6FpO0RsV7S9ux5IzXhS8xfws3QGSr8n6EOcoWM7bWS/mdJH+uYvFnSddnj6yS9aaotS8QsfJFn4T02Df9nqIu8lcwfSvo/JR3pmHZKRCxIUvZzzXSbNpm5AwsTr6NpX+SmvR8A6RsaMrZ/QdLeiPjGOBuwfZntHbZ37Nu/f5xVVGJWdsiz8j4BVCNPJfNqSb9oe5ekGyS93vafSHrY9rwkZT/39lo4IrZGxMaI2Lh61aopNbtYTd7xdr63Jr9PAGkYGjIRcWVErI2IdZIukvQ3EfGrkrZJ2pLNtkXSZwprZWd7Ej1Bs044IglAWSY5T+YqSefbvl/S+dnz2mPnC6AItjfZvtf2TttHHY3rlg9lr99p++xhy/Y7lcT2+ba/Yfvb2c/XdyzzK9n677L9/q42vNX23dlrfzaN9z1SyETElyLiF7LH+yPivIhYn/18dBoNAoCmsT0n6WpJF0jaIOli2xu6ZrtA0vrs32WSPpJj2X6nkuyT9M8j4h+r1dP0iWxdqyR9QNJ5EfEKSafYPi97bb2kKyW9OnvtN6bx3mt1xn/RqGIAFOQcSTsj4oGIeEqt8e3NXfNslnR9tHxV0gnZePegZXueShIRt0fEQ9n0uyQ9z/axkl4s6b6IeCR77QuS3pw9/t8lXR0Rj2Xr6DnOPqrahUyKNy4DgCFOlfSDjud7sml55hm0bJ5TSd4s6faIeFLSTkkvt73O9lK1Qum0bL6XSXqZ7b+3/VXbm0Z7i7019lL/iyvnq24CgNmy2vaOjudbI2Jr9tg95o+u5/3mybNsT7ZfIek/SHqjJEXEY7bfKemTap33+N/Vqm6kVh6sl/Q6SWsl/Tfbr4yIx/Nsq5+kQ6b7pmVFVjF0lQHN42XHatn8urI2ty8iNvZ5bY+erRik1k78oZzzHDNg2Ydtz0fEQvepJNmVWj4t6W0R8b329Ij4rKTPZvNcJmmxY/tfjYjDkr5v+161Que2ge96iNp1lwFADd0mab3tM2wfo9bpINu65tkm6W3ZUWbnSvpR1gU2aNmep5LYPkHSX0m6MiL+vnMjttdkP0+U9K/17OXC/lLSz2WvrVar++yBSd94EpXM4sr5oZeBGaWKoasMQEoi4mnbl0u6RdKcpGsj4i7b78hev0bSzZIuVGvc5JCktw9aNlv1VZJutH2ppN2S3pJNv1zSSyX9ju3fyaa9MRvM/6Dtn8mmvTci7sse3yLpjbbvVqu6+T8iYuLLtCQRMoN0d5kVga4yAEWLiJvVCpLOadd0PA5J78q7bDZ9v6Tzekz/PUm/12ddF/eZHpL+TfZvamrRXcYRZQBQT7UImVHQVQYA6WhcyIyKrjIAKE7yYzKjoIpBmVYcOfSc5/zBAhwtmZDJc4QZUJXuQBk0D2EDPCuZkKkCOwP0kydUhi3L5wtoUMjQVYZJTRIs/dZH0GDWNWbgn642TGragUDAAA0KGWAaCAZgupIJGSoRpIKgAaYnmZCZBoIK03JwyfKJwoagAloaFTLAtBEWwGQIGWCIUYOGYAKelXzIxN7dI12FmS4zFIHgAMaTRMgQDKiDPEFDGAHPlUTI5FH0PWWAPAgRYDRJh8y4wUJlhCJNeuQZMEuSDpluVDNISXfQEDzA0WoVMqOoUzWz4sih5/xrmqa+L4lgAYap/AKZ/cKgX9USe3c36nbMvXa+vaY1YWfW+b6a8H7amvRegGmrPGSQT9N20FwOH5gNje0uAwBUr/KQ6XcfmH5dYk3qKhtH3f7yb+pYDIB86C5DpbixF4rkpcfM/B+mVau8khkFHxYAqJekQ4ZQAYB6SzpkOhE4AFA/SYRMv8H/pmv6oHjT3x+A4ZIf+KeCAYD6SqKSKcKsVkd1RMUDNFdjQwYAUD1CBoWgOgEgJRQydG8Nx0mLAOommZABADQPIQMAKAwhAwAoTFIhM61xmTzrYXwjHfxfAM2VVMhMAwcQpCFvcBAwQLM1KmQImHohYIDmSy5kygiKOu7cmtbmOr4fAKNLLmTGRRVTHwQMMDsaETIETLq6A4WAAWZLkiEzSmiMGjB5dnL7n1gcaZ3Ih4ABZk/yl/ofhAqmHggXYHYlWckAAJoh2ZAZVqUUXcWk1GVGJQCgrpINmUHGDRjGYwCgXLULGcZhAKA+kg6Z7kApO2CoagBgMkmHzDQxrgEA5av1IcyzgHAExhdzy+hir1jSlczcgYWqmwAAmEDSIVM2xmAAYLpqFTJUNgBQL8mGDIECAPWXbMikgi40ABhf7UJmliocjiwDUHdJhswsBQkANFmSITMMIQQA9VDLkAEA1ENyIUOVAgDNkVzI5DVKGHGJfwCoRm1DBgCQPkImURy+DKAJkgoZxmMAoFmSChkuyQ0AzZJUyAAAmqW2IdPkqofxGABNkVzINDk8AGDWJBcyVVl13FzVTQCAxqllyDSl2unVLUZXGYAmSTJkUgoRKhwAGN/QkLH9PNtft/0t23fZ/t1s+km2b7V9f/bzxOKbm1YAAQAGy1PJPCnp9RHxM5LOlLTJ9rmSrpC0PSLWS9qePccE6CoDmsv2Jtv32t5p+6j9pVs+lL1+p+2zhy076I9921dm899r++c7pv+s7W9nr33ItrPpx9r+ZDb9a7bXTeN9Dw2ZaDmYPV2W/QtJmyVdl02/TtKbptGgNioWAE1he07S1ZIukLRB0sW2N3TNdoGk9dm/yyR9JMeyPf/Yz16/SNIrJG2S9OFsPcrWe1nHtjZl0y+V9FhEvFTS/y3pP0zjvecak7E9Z/sOSXsl3RoRX5N0SkQsSFL2c800GjTIOMFDdQAgAedI2hkRD0TEU5JuUOsP9U6bJV2f/WH/VUkn2J4fsmy/P/Y3S7ohIp6MiO9L2inpnGx9KyPiKxERkq7vWqa9rpskndeuciaRK2QiYjEizpS0NmvoK/NuwPZltnfY3rFv//4xm9lchCAwE06V9IOO53uyaXnmGbRsvz/2B61rT591PbNMRDwt6UeSVuV6dwMsHWXmiHjc9pfUKq8etj0fEQtZOu7ts8xWSVsl6WfPOjNG2d7iyvlnLprZ9O4zwgaYvkUvKfO7tdr2jo7nW7P9nyT1qgi694f95smzbLdx1jXOdobKc3TZybZPyB4fJ+kNkr4raZukLdlsWyR9ZtLGpIbDlwGMYF9EbOz4t7XjtT2STut4vlbSQ13L95tn0LIPZ3/kq+uP/UHrWttnXc8sY3uppBdIenTQG84jT3fZvKQv2r5T0m1qjcn8V0lXSTrf9v2Szs+eYwxUMUDj3SZpve0zbB+j1qD8tq55tkl6W3aU2bmSfpR1gQ1att8f+9skXZQdMXaGWgP8X8/W92Pb52bjLW/rWqa9rv9F0t9k4zYTGdpdFhF3Sjqrx/T9ks6btAHD1KmbbMWRQwQGgKNExNO2L5d0i6Q5SddGxF2235G9fo2kmyVdqNYg/SFJbx+0bLbqqyTdaPtSSbslvSVb5i7bN0q6W9LTkt4VEe17zL9T0h9LOk7SX2f/JOnjkj5he6daFcxF03jvnkJQ5Xbm2WfH17/4+dK2J41eJex/YvGZx6N2l604cmisbQI42qqVx38jIjZOso4zzz47/uZv/25aTRpoGu1totIvK9PUHXA7YLofA8Asq+TaZQeXLG9s2LQRNABQ8QUyUw6aaRxZRtAAmHWVX4U55aABAEym8pCR0gqaaZ8bQzUDYJYlETJSWkEzbQQNgFmVTMhIBA0ANM1I1y4rQztoJt0pTxJY454fAwB4ruRCpu3gkuW5d951qYC4IgCAWZNsyEj9g6bOO2qCBsAsSTpkpHoHSj8EDYBZkdTA/yxhHAfALCBkAACFIWQAAIUhZAAAhSFkAACFIWQAAIUhZCbEUWIA0B8hAwAoDCEDAChM8iGz/4lF7X9isepmAADGkHTIdIYLQQMA9ZNsyPQKFYIGAOolyQtkDgqT/U8sTv0WyQCaafEIf5xWLdlKZhA+NABQD8mFTN4AIWgAIH1JhUy/4PjhwadKbgkAYBqSCZlxAoZqBgDSlkTI5AmYfmFT16DhzpgAZkESIdNLr1BpWtAAQNMlGzIAgPprdMhwhWQAqFYjQqbXyZntgCFoAKA6lYcM4ykA0FyVh0wvgw5bHuecGaoZAKhGkiEDAGiGmQkZqhkAKF8jQ4ZAAYA0JBcyRV6njPABgHIlFzLddj3+hHY9/kTVzQAAjCH5kJk2qhkAKE8tQ4ZL/wNAPVQaMsNOxCyqm4xqBgDKkVQlM6hCyRs4BAgApCOpkAEANMvMhgwVDwAUL9mQafphy4QcgFmQbMhMw9yBhULXT1AAwGDJhMy0D0suOmAAAMMlEzJ5jNOFRtgAQHWWVt2AXjrDZOcjB/XSk1cUsp0VRw7p4JLlhawbQPUOLx7h5O2K1aKS2fnIwaOm8cEBgPTVImQmRZcZAFSjspAZdkkZAED9JVHJdHZ9dY/HjIJDigEgLUmETB7twMlzhFmv7rF+XWYEEwAUpzYhAwCon2RDZtSusjw4AAAAypVUyPTqCvvuwoEKWgIAmIakQqYfggYA6qnykOl1UmW/rrJRBv8lKfbuPmpary4zBv8BoBiVh8wgk1QwvQIGAFCupEMGAFBvyYTMOFdYznP9MioaAKhOMiEzDIP/AFA/yYVMEefHpIhbDACYBcmFTFu7ctnzw9FCJ88Jl5yUCQDlSDZkOnUGzTiVDuMyAFCNWoQMAKCeKgmZptxLhnEVABis0kqmfQhy+/DldlfYuOMxAIC01Kq7rPMw5nHOqwEAlKtWITMJBv8BoHzJh8zje38iabSuszyBwmHMAFJh+yTbt9q+P/t5Yp/5Ntm+1/ZO21fkWd72ldn899r++Wzactt/Zfu7tu+yfVXH/C+yvd32nba/ZHttx2un2/687Xts32173bD3lmzIdIZKO2hGdXhh15RaAwCFukLS9ohYL2l79vw5bM9JulrSBZI2SLrY9oZBy2evXyTpFZI2Sfpwth5J+oOIeLmksyS92vYF7emSro+IV0l6r6Tf72jG9ZI+EBE/LekcSXuHvbHkQibP5WPynCtDwACokc2SrsseXyfpTT3mOUfSzoh4ICKeknRDttyg5TdLuiEinoyI70vaKemciDgUEV+UpGxd35TUrlg2qBVUkvTF9jaywFoaEbdmyx2MiKH3SUkmZMq4nAzjMgASdUpELEhS9nNNj3lOlfSDjud7smmDlh+0jCTJ9gmS/rmeDZZvSXpz9viXJD3f9ipJL5P0uO2/sH277Q90VEV9LR02Q5V6dZN9d+GAXj6/soLWAKibJxePlHkk6mrbOzqeb42Ire0ntr8g6YU9lnt3zvW7x7SYZBnbSyX9uaQPRcQD2eR/K+k/275E0pclPSjpabXy4p+p1b22W9InJV0i6eODGjA0ZGyfplY/3AslHVHrF/dB2ydlG1knaZekt0bEY8PWV7bDC7u0bH5d1c0A0Hz7ImJjvxcj4g39XrP9sO35iFiwPa/eYx17JJ3W8XytpIeyx/2WH7SMJG2VdH9E/GFHOx+S9MtZu1ZIenNE/Mj2Hkm3t8PI9l9KOldDQiZPd9nTkn4rG+g5V9K7sr65oQNVefT6K6PqkzCbckUCALWxTdKW7PEWSZ/pMc9tktbbPsP2MWoN6G8bsvw2SRfZPtb2GZLWS/q6JNn+PUkvkPQbnRuxvdp2OxuulHRtx/ZPtH1y9vz1ku4e9saGhkxELETEN7PHP5Z0j1p9enkGqqaq6vABgIJcJel82/dLOj97Lts/ZftmSYqIpyVdLukWtfbDN0bEXYOWz16/Ua0w+Jykd0XEYnZY8rvVGuT/pu07bP9v2bpeJ+le2/dJOkXS+7J1LarVlbbd9rfV6or76LA3NtKYTHZM9FmSvqaugSbbvQaqAABDRMR+Sef1mP6QpAs7nt8s6ea8y2evvU9ZUHRM26Pe4zWKiJsk3dTntVslvarf++gl99FlWd/cpyT9RkTkvk2l7cts77C9Y/++faO0TZL0+AInTQJAXeUKGdvL1AqYP42Iv8gmP5wNMGnAQJUiYmtEbIyIjatWrx67oeOekAkAqM7QkLFttY4euCci/lPHS3kGqkaS50TMcXBiJgBUI08l82pJ/1LS67PBoTtsX6g+A01FmGaXGSdkAkB5hg78R8Tfqc8AkfoMNE3DqN1jK44MvbrBUeYOLGhx5fzIywEA8knmsjKp4VwZAJgcIQMAKExlIdO+9fI0dd8j5tCuXVPfBgAgvyQqmc4rMD/nPjKcIwMAtZZEyExT7N3d85BlDmMGgPI1LmTa6CoDgOo1NmQAANVLPmQOPHhf1U0AAIwp+ZABANRX40NmkrEZTsgEgMk0PmR64fplAFCOmQwZAEA5CBkAQGGSDBluUAYAzZBkyAAAmoGQAQAUhpABABRm6J0xm+Twwi4tm19XdTMAlOTJp4885yrvKB+VDACgMIQMAKAwhMwQXFoGAMY3EyEzyvXLCBUAmJ6ZCBkAQDUIGQBAYQgZAEBhCBkAQGEIGQBAYZIMmRPWHF91EwAAU5BkyAAAmiH5kFl56suqbgIAYEzJh0wTHVyyvOomAEApCBkAQGEIGQBAYQgZAEBhCJkcuGgmAIyHkAEAFIaQAQAUhpABABQmuZBZ+8IVzzw+YX6+wpYAACaVRMi89OQVw2cCANROEiEDAGimykLmhSuOKW1by9etK21bAIBnzVQls2x+XdVNAICZkmzIFHlPGa85vbB1AwCetbTqBuRxwvy8Tlhz/HOOPJuGxZXlH73GFZiB8vx/hxf13YUDVTdjpiVVybx8fuVU5gEApCGpkAEANEvyIdPrhMxpn1dDFxYAFCP5kAEA1FflIbPuhOOOmjbNAX7OkQGA6lQeMuPqFU6po1sOwKxJOmSmea5MnhMxuTkZAExX0iEDAKi3xoZMv7GYcc/2p8oBgNElEzJc7h8AmieZkAEANE9yIdPrsjGjHgDQHuTn8GUAqFZyIVOWsi+OyeHLAGZRsiHTPiGzXcX0OkGzzBufAQBGl2zI9JLnCsy9jh4r8mZlVCgA0F+tQmZS3KwMAMpVaci0u7val4hpH8bcrlh6dZHV8QrMVDsAZlUlIbPquLlC1lvFnS4BAP3VortskqsyFzkeAwAYLNmQKfI2y1Q8AFCO5EKme8ylXcWMGjrdFUxVg/6MxwCYZcmFTB51vJcMAMyiZEKmiOBgPAZAHdg+yfattu/Pfp7YZ75Ntu+1vdP2FXmWt31lNv+9tn++Y/qXsml3ZP/WZNNfZHu77TuzedZm08+0/RXbd2Wv/Uqe95ZMyPQy7rgM58MAqJkrJG2PiPWStmfPn8P2nKSrJV0gaYOki21vGLR89vpFkl4haZOkD2frafsXEXFm9m9vNu0PJF0fEa+S9F5Jv59NPyTpbRHRXtcf2j5h2BtLOmS6jXuOTNGh02/chfEYADltlnRd9vg6SW/qMc85knZGxAMR8ZSkG7LlBi2/WdINEfFkRHxf0s5sPYNsUCuoJOmL7W1ExH0RcX/2+CFJeyWdPOyNVR4yva4/1h0m0zzSjCPLACTolIhYkKTs55oe85wq6Qcdz/dk0wYtP2gZSfqjrKvsd2w7m/YtSW/OHv+SpOfbXtXZENvnSDpG0veGvbGlw2YAAOSy2vaOjudbI2Jr+4ntL0h6YY/l3p1z/e4xLSZY5l9ExIO2ny/pU5L+paTrJf1bSf/Z9iWSvizpQUlPP7NCe17SJyRtiYgjwxqdfMgUeb4M3VlAsz11+Ij2/PBgWZvbFxEb+70YEW/o95rth23PR8RCthPf22O2PZJO63i+VtJD2eN+y/ddJiIezH7+2PafqdWNdn3WFfbLWbtWSHpzRPwoe75S0l9J+u2I+Gq/99Op8u6yTk06NJkAAzCCbZK2ZI+3SPpMj3luk7Te9hm2j1FrQH/bkOW3SbrI9rG2z5C0XtLXbS+1vVqSbC+T9AuSvpM9X227nQ1XSro2m36MpE+rFUT/Je8bSypkOk3rQpgcaQagBq6SdL7t+yWdnz2X7Z+yfbMkRcTTki6XdIukeyTdGBF3DVo+e/1GSXdL+pykd0XEoqRjJd1i+05Jd6jVJfbRbF2vk3Sv7fsknSLpfdn0t0p6jaRLOg57PnPYG6usu2zVcXPa/8SipNbg/w8PPjVw/nbotKud9gED7YttHlyyXCuOHCqquYVd1BMAImK/pPN6TH9I0oUdz2+WdHPe5bPX3qdng6I97SeSfrbP/DdJuqnH9D+R9CeD3kcvyVUynV1m41Qz7aPHelUwHFkGAOVKJmTy3Eq5u4oZhuuVAUC1Kg2Zfl1Q3dXMtG9UNmp76CoDgPEMDRnb19rea/s7HdNyXWdnVJ0VSp4jzUbZ+RfdVUb1AgBHy1PJ/LFa16npNPQ6O9M2alfZtIxaxRA2APCsoSETEV+W9GjX5DzX2cmleyc+ajXTS56qhTAAgOKNOyaT5zo7UzNpFTNuVxlVDABMpvCBf9uX2d5he8f+fftyLZOnmilqMJ5BfgCYnnFD5uHs+jjti6X1us6OJCkitkbExojYuGr16p7z9Nqxlz320o2wAYDJjRsyea6zM1Vld5UBACaX5xDmP5f0FUn/yPYe25eqz3Vypq3qagYAMJmh1y6LiIv7vNTzOjnj6ryWWbdxwmaSCoauMgCYjmQuK9NPr4AhBACgHpIPGQBAfSUVMlQoANAsSYVMWThpEgDKUdlNy/ppVzO9DgIYVOkUecMyAMB4kguZtu6w6Rcwo4YLVQwAlCfZkGmbZvVCwABAuZIPmV7G6RojYACgfLUKmXHHXQgYAKhGbY4uI2AAoH6SD5kVRw6VfuRYFdvj6DgATZR0d9mkO95xqpj2NlccOVRKFdT5HtuPqb6A6Vg8fESP7/1J1c2YacmGTBV/2Xdvs+ig6fceiwqbftsj1AAUJbmQSa3bqKigyfM+pxE2o2ynE8EDYBqSCpleO7u5AwuSRr90/6g7yUE742lXFqMG6ajbn0ZQd69j2Lan9ccB4QY0SzIhMyhg2o+Lustl3h3kNKqaSXbGnct2t6PoCrCsCrOMsBq2DYIOmJ7KQ6bfF74zYDqn5QmaUXYS41QV4+6EprmjTq1bMTXTCvNeCCEgv0pDZlj1kqpxgoZQaI4i/i8JLjRVZSGTN2Bi725Jktec/sw80+o2m/Sv3TLHSNBsVE9oqkpOxhw1YLofD6p2Jt3xzx1YyF1NjXvkFjCq9gm73f+A1JVayczFkbECpue6CjwQoL1+afhRbYOO/GIngKJx+DlSV/llZUYJmO7p447fjHKwwbgIGFSFigcpqTxkUjdO1xlfaqSGwEFVKg2ZPFXM4YVdOrywq+/rAEZD2KBMja5k2rduntQo3Wh8eVEXhA3KkFTI9Kpiej0uwjTGY/jCoo4IGxQpqZDJK+/hzNNUh5NEgUkQNihCsiHTrlwO7dqlQ7t2VdoWYJYQNpimykKmuzLIM6BfdJcZgGcRNpiGJCuZziomJXSZYRYRNphEkiHTS3fg5D2UeVpHmAGzjrDBOJIImc7A6FXF/Oh7Dx71+jCjhMsoFQrVDGYdQYNRJBEyg3QGTL/us2nv+DnhExiMqgZ5JRUyKYzFtAOGoAGGI2gwTCUh01l5DNqZd1YxvbrMqgoCusyAZ1HVYJBkKpm8Vcw4Vc64g/9UM0B+BA16qTRk8lQxj+18RI/tfOSo16s+Z4ZqBjgaVQ26JVHJ5K1iOrvM2iatNggLYPoIGrSVemfMvDqrmLbHdj6iE196sqRWGC1ft27q2+0VWLF3t7zm9J7zF313TqDOBt01tiyLh5/S4wv8IVmlJCqZaaAiAYD0NCZkAADpqXXITGPwv7O7q1e3WL+uMgDAcKWHzKTdWr0G/wEAaaqskmkPsg+rRh7d+age3fnoVLc9rYFIBv0BYLBkusvahy/3OrKsKnSVAcBkkgkZAEDzEDIAgMLUKmQ6u9A6rw5QxDXG6CoDgMnVKmTymPRyFoQLAExP40JmGggaAJiOpENm2ocu98OhyABQjCRDptfhy2UFDgBgepIMmWHGOet/3BuXAQDGl1TIcMkYAGiWpEKmnwf39j9irOo7ZAIA+ksiZAbdEXNQwAAA0pZEyIwihWuaAQDyqV3ItG/BPKpVx80953n3lZg5jBkApi+JkFm+bl3f105ds1ynrqnuHuH9EEoAMFwSIdP2gpecWnUTciFgACCfpEKmbViXWGcYLZtf13Oe7u6xaSFgACC/JEOm7aSXntTzcVEGBcjiynkCBkAhbJ9k+1bb92c/T+wz3ybb99reafuKPMvbvjKb/17bP98x/WLb37Z9p+3P2V6dTT/d9hdt3569dmHHMu+3fZfte2x/yLaHvbekQ2ZaJq1qCBcABbtC0vaIWC9pe/b8OWzPSbpa0gWSNki62PaGQctnr18k6RWSNkn6sO0520slfVDSz0XEqyTdKenybF2/LenGiDgrW/bD2br+R0mvlvQqSa+U9D9Ieu2wN1arkBn3yLJJEDAASrBZ0nXZ4+skvanHPOdI2hkRD0TEU5JuyJYbtPxmSTdExJMR8X1JO7P1OPt3fFaNrJT0ULZMZM8l6QVd058n6RhJx0paJunhYW+sViFTNgIGQElOiYgFScp+rukxz6mSftDxfE82bdDyPZeJiMOS3inp22qFyAZJH8/meY+kX7W9R9LNkn4tW+9XJH1R0kL275aIuGfYG1s6bIaieM3pir27tWx+nQ4v7NLydet0aNcuveAlp+pH33tQJ7705IEnXg467Dmvg0uWH3WTs8WV85o7sEDAAA2wePhJHXjwvrI2t9r2jo7nWyNia/uJ7S9IemGP5d6dc/29xj9inGVsL1MrZM6S9ICk/0fSlZJ+T9LFkv44Iv6j7X8q6RO2XynpxZJ+WtLabD232n5NRHx5UANKD5n2TjyvMgb8uxEwAMawLyI29nsxIt7Q7zXbD9uej4gF2/OS9vaYbY+k0zqer9WzXVn9lu+3zJlZm76Xbf9GPTsOdKla4zeKiK/Yfp6k1ZJ+SdJXI+JgtsxfSzpX0sCQqV13WV3OpQGAEWyTtCV7vEXSZ3rMc5uk9bbPsH2MWoPy24Ysv03SRbaPtX2GpPWSvi7pQUkbbLcHus+X1O762i3pPEmy/dNqjcM8kk1/re2lWSX02o5l+qqsu6wsRZ0vAwBTdJWkG21fqtbO/C2SZPunJH0sIi6MiKdtXy7pFklzkq6NiLsGLR8Rd2VVyt2Snpb0rohYlPSQ7d+V9GXbhyX9g6RLsnX9lqSP2v5NtbrjLomIsH2TpNerNY4Tkj4XEZ8d9saSDpnOcZkqjiwDgDJExH5l1UPX9IckXdjx/Ga1BuNzLZ+99j5J7+sx/RpJ1/SYfrdahyp3T1+U9K8GvY9ekuwuo0sMAJohqZCZxhFjAIB0JBUyVei+5D8AYHoaETJec3rVTQAA9JDEwH/7hMxeOgf822M17W61fldgbuPIMgCoViWVTPtkx14VSDtAGPwH6o/uaJQaMos+enPtoOlVlXQGTb8qpr18O7j4UANAOkqvZNoh0Hnplu6g6TzK7AUvOZWAAWqI7yOkirrLegVNW6+gKas9AKaD7xTaJgqZfndpG0We8ZnOx1QxQNr4LqLT2CEz5C5tQ3V+ELuDpnN8Zvm6dQQMUBN8F9Ftkkpm0F3acun1gRx0IAAAoF4mCZlBd2nLLc+BAJ2Pi6pi+AsMmAzfIfQyScjkukub7cts77C9Y/++fT1XNOxAALrJgLTxXUQ/k4TMoLu0PSMitkbExojYuGr16qErHXQgAID0EDAYZJKQGXSXtpENOhCg8zFVDJAOvocYZuxrlw25S9tYDi5ZrhVHDj1nWmfQEDAAUC8TXSCz313aJtEOmsWV85o7sDDNVefeNoDh+GMPeSR5qf9eBwJQxQDp4HuIvJK41P8gnUHDBxuoHt9DjCLJSkbigwykiO8lRpVsyEjP/UCX9eHmSwQA0+OIo86fLG5j9o8l3VvaBgdbLan32aHlS6ktUlrtoS29pdQWqZj2vCgiTh4+W3+2P6dW28qwLyI2lbSt2ig7ZHZExMbSNjgAbekvpfbQlt5SaouUXnuQjqS7ywAA9UbIAAAKU3bIbC15e4PQlv5Sag9t6S2ltkjptQeJKHVMBgAwW+guAwAUppSQsb3J9r22d9q+ooxtDmnPLtvftn2H7R0lb/ta23ttf6dj2km2b7V9f/bzxArb8h7bD2a/mztsX1hSW06z/UXb99i+y/avZ9Or+t30a0/pvx/bz7P9ddvfytryu9n00n83A9pSyecG6Su8u8z2nKT7JJ2v1j1obpN0cUTcXeiGB7dpl6SNEVH6eQa2XyPpoKTrI+KV2bT3S3o0Iq7KQvjEiPh3FbXlPZIORsQfFL39rrbMS5qPiG/afr6kb0h6k6RLVM3vpl973qqSfz+2Len4iDhoe5mkv5P065J+WSX/bga0ZZMq+NwgfWVUMudI2hkRD0TEU5JukLS5hO0mKSK+LOnRrsmbJV2XPb5OrZ1ZVW2pREQsRMQ3s8c/lnSPWrfzrup30689pYuWg9nTZdm/UAW/mwFtAXoqI2ROlfSDjud7VNGXtUNI+rztb9i+rOK2SNIpEbEgtXZuktZU3J7Lbd+ZdaeV0j3VyfY6SWdJ+poS+N10tUeq4Pdje872HZL2Sro1Iir73fRpi1Tx5wZpKiNk3GNa1X/5vDoizpZ0gaR3Zd1GaPmIpJdIOlPSgqT/WObGba+Q9ClJvxERB8rcds72VPL7iYjFiDhTrducn2P7lWVsd4S2VPq5QbrKCJk9kk7reL5W0kMlbLeviHgo+7lX0qfV6tKr0sPZGEB7LGBvVQ2JiIezncgRSR9Vib+brI//U5L+NCL+Iptc2e+mV3uq/P1k239c0pfUGgOp9HPT2Zaqfy9IVxkhc5uk9bbPsH2MpIskbSthuz3ZPj4byJXt4yW9UdJ3Bi9VuG2StmSPt0j6TFUNae+0Mr+kkn432YDyxyXdExH/qeOlSn43/dpTxe/H9sm2T8geHyfpDZK+qwp+N/3aUtXnBukr5WTM7HDGP5Q0J+naiHhf4Rvt35YXq1W9SK2btv1Zme2x/eeSXqfWlWEflvR/SfpLSTdKOl3SbklviYjCB+T7tOV1anV5hKRdkv5Vu9+/4Lb8T5L+m6RvSzqSTf73ao2DVPG76deei1Xy78f2q9Qa2J9T6w/DGyPivbZXqeTfzYC2fEIVfG6QPs74BwAUhjP+AQCFIWQAAIUhZAAAhSFkAACFIWQAAIUhZAAAhSFkAACFIWQAAIX5/wG8kUTu7iPQfQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "zeta_bt = (zeta * ds.e3f).sum(dim='z_c')\n", "plt.contourf(\n", " domcfg.glamf,\n", " domcfg.gphif,\n", " zeta_bt.isel({'t':0}),\n", " levels=np.linspace(-abs(zeta_bt).max(), abs(zeta_bt).max(), 21),\n", " cmap='RdBu_r'\n", ")\n", "plt.colorbar()" ] }, { "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", "$$ \\Psi(x,y) = \\int_0^x \\int_{bottom}^{surface} v_{bt}(x,y) \\, \\text{d}z \\, \\text{d}x $$\n", "\n", "We calculate this by integrating $v_{bt}$ along the X axis using the grid object's `cumint` method:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Coordinates:\n", " * t (t) object 1050-07-01 00:00:00\n", " * y_f (y_f) float64 0.5 1.5 2.5 3.5 4.5 5.5 ... 35.5 36.5 37.5 38.5 39.5\n", " * x_f (x_f) float64 0.5 1.5 2.5 3.5 4.5 5.5 ... 15.5 16.5 17.5 18.5 19.5\n" ] }, { "data": { "text/plain": [ "(0.0, 60.0)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAJDCAYAAAALuGNCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6fUlEQVR4nO3de7RdZ1nv8d+TpE1pk7Zp06YhbQmVeqPaFmK94FGkgBXRolCEM/BEZZyqQxjgZUgROXjBM+qNQ496POZQsCgClYvtEAdYquDlHCspVKAERGNN06a5lJRkt5K2yXP+WHOlM2uvy7zP933n9zPGHnuvudac883O2u9vPe87L+buAgBgRd8NAACEgUAAAEgiEAAAGQIBACCJQAAAZAgEAICkgoFgZmea2fvM7PNmtsPMvtXMzjKz28zsi9n3dW03FgDQnqIVwg2SPuzuXyvpUkk7JF0n6XZ3v1jS7dljAECkbNGJaWZ2uqR/knSR515sZl+Q9Gx332NmGyV9zN2/ptXWAgBaU6RCuEjSfknvMLNPmdnbzOw0SRvcfY8kZd/PbbGdAICWrSr4mmdIerW732FmN6jE8JCZXSvpWkl60qmnPfOip11cqaEAhuXuT991wN3PqbONi1ac6o/oaFNNmmuvP/oRd7+qk521pEgg7Ja0293vyB6/T6NA2GtmG3NDRvumrezu2yRtk6RLLr3cb/7IxxtoNoDUPX3jGf9edxuP6Ki2rtrURHMW+o3H/m19Jztq0cJAcPcHzOxeM/sad/+CpCslfS772irp+uz7La22FEB01p2ysu8moIQiFYIkvVrSu8zsZEk7Jf2oRvMPN5vZKyXtknRNO00EAHShUCC4+12Stkx56spGWwMA6A1nKgMAJBEIAIAMgQAAkEQgAAAyRY8yAtCxB5aOzH3+vDWrO2oJhoJACNiiDgHDln9/EA5oAoFQE502QjD5PiQgUAWBsAAdPmL0wNKR2qFAyAwPgZBD54+U1AkF/haGafCBwBsfAEYGFwgEAIakSpXA38hwDSIQeIMDxfC3MmzJBgJvbKBZTUxUI2xJBQIhACxXtCPn7wdJBAJvZKAe/oYgRRwIIb+B7zsUbtvasOl0hhFCx3BPv8zsAknvlHSepGOStrn7Dbnnf07Sb0o6x90P9NPKyAIhlBAYWoe/SF+/D4KoGaH8XSXucUk/6+6fNLO1ku40s9vc/XNZWDxPoztP9iqaQOjzTUsAhGnW/wtBUQxB0B133yNpT/bzYTPbIWmTRvem/x+Sfl4B3Jc+ikDo+o1LAMSt7P/fUAKk7t8RQ07NMLPNki6XdIeZfb+k+9z9n8ys34YpkkDoEmEwPFQa6Vq9wrT51JO72dmXtd7MtueWbHP3bfmXmNkaSe+X9FqNhpHeIOn53TRwseADoavqgCDApPx7gnBAAQfcfcusJ83sJI3C4F3u/gEz+wZJT5U0rg7Ol/RJM7vC3R/opMUTgg+ELhAGWGTyPUJAoAwb9fg3Strh7m+RJHf/jKRzc6+5R9IWjjKaoe3qgCBAVVQPKOlZkn5Y0mfM7K5s2S+4+1/016Tlgg6ENhEGaMr4vUQwYBZ3/ztJc2eN3X1zN62ZLdhAaKs6IAjQFoIBsQs2EJpGEKArBANitaLvBnSBMEAf7jt0JJn3HucgDEOQgdDkcFEqf5CIF+9BxCLIQGgKf4gIRUrVAtKVbCDwx4cQEQwIWbKBAIQspmBg/mA4kgyEWP7QgJiCAelL7rDTUP64dh18pNHtXbju1Ea3h7BwqCpCEFwg1DnCqM8waDoAmtg+IRKf+w4dIRTQm+ACoao+wqDtEKhrUfsIjDCFVC0wfzAsyQRCl0IPgqLm/TsIi/6FFAwYhiQCoavqIJUgKIKwCAfBgK5EHwhdhMGQgqCIyd8HAdEN5hfQtugDoW2EwWIERHcIBbQp6kBoszogCKojINpFKKAtQQVCV/dPXoQwaFb+90k4NINQQBuCCoQy2qoOCIN2UT0A4Yo2EJAGAqI6qgQ0LcpAoDpIF8NL5RAK8528YoUuXHtyNzv7cje7aVOUgdAGwiA8hAPQregCoY3qgDAIH+EwG1UCmpLk5a+Rtl0HHyHEgRZEFQhUB8gjGJ4QymXfEbeoAqFpdCZp4P8RaMagAwHpIBTaEcrJouhGNIHQdElMB5IehpCAeqIJBKCooYYC8wioK5hA6LI0HWqHMST8HwPlBRMI8zT5yYeOYjj4vwbKiSIQgKoIBaC46M5UriO0zmHn/odLr3PROae10JK07Tr4CGc3AwUMKhD6ViUAim6DoJhvKKHAZSxQx2ACoY/qoIkAqLovAmK5oYRC0x5YOqLz1hAyQzCYQOhKlyEwz7gdBAOAogYxqdxVdRBKGOTt3P9wkO0CEJ7kA2HIYZBHMABYJPhAiOHsy5g62pjaCqBbwQdCHV1UBzF2sDG2GUD7kg6EtsXcscbcdnSPq54OA4FQUQodKvMKaYphmBVhSjYQ2hwuSq0TTe3fM09oZ6sDIUk2ENqSaueZ6r8LQHGcmFbQEDrMnfsf5kQ2JGXlySt15lPO6GZnu7vZTZuSrBCaHhYYQhiMpf5v5dIVwGxJBgIwZFzcDlURCAuk/ol5SKgOgPmCCIQmj3FucriIMEBsqA5QRxCBMEufx1MTBmmhOqiPk9PSF3Qg9IUwQIyoDlBXUoHASUf1pXjYKdUBUExSgdAEqgMAQ0UgIGlDqQ4YLkITkgmEJoaLqA4ADFkygQBMojoAyiEQMlQHAIYuiUDg6KJmpHSEEdUBUF4SgVAX1QFhECPCAE0rdPlrM7tH0mFJRyU97u5bzOwsSe+VtFnSPZJe6u4H22km2kQYxKePMDhvDQGUujIVwne5+2XuviV7fJ2k2939Ykm3Z487V3e4iOogHYQBUE+dIaOrJd2U/XyTpBfVbk0O94XtRirVAWEA1Fc0EFzSX5rZnWZ2bbZsg7vvkaTs+7ltNLBNQ68OUgmDoSAM0LaigfAsd3+GpO+R9FNm9h1Fd2Bm15rZdjPb/qUHH6zUyFk4uqi6lMJgCNUBYRA/M3u7me0zs8/mll1mZv9gZndl/eQVfbaxUCC4+/3Z932SPijpCkl7zWyjJGXf981Yd5u7b3H3LWedfXYzrW7A0KuDVBAGiMgfSrpqYtlvSPpld79M0n/LHvdmYSCY2Wlmtnb8s6TnS/qspFslbc1etlXSLW01chqqg+pSqQ4IA8TE3f9G0pcmF0s6Pfv5DEn3d9qoCUUOO90g6YNmNn79n7j7h83sE5JuNrNXStol6Zr2mtmsIVcHhAGq4JDT1rxW0kfM7Lc0+oD+bX02ZmEguPtOSZdOWf6gpCvbaBQAqoMmrFq9UusuOrObnf291pvZ9tySbe6+bcFaPynpp939/Wb2Ukk3Snpua21coNCJaaGpM1xEdRC/IVQHhEGUDuTO0ypqq6TXZD//qaS3Ndukcrh0xUAQBvEILQwYLmrV/ZK+M/v5OZK+2GNb4qwQqhpqdZBCGAwhCKTwwgDNMbN3S3q2pPVmtlvSmyT9V0k3mNkqSV+RdO3sLbSv90B4YKncGckcXVRcCkEgDScMQkR10Bx3f/mMp57ZaUPmCHLIiMtW1EcYxIfqAH0LMhDaMKThIsIgPiGGAdXB8EQVCAwXLUYYxIcwQCiiCoSqhlIdEAbxCTEMMFzRBALVwXyEAZpCdTBc0QQCZiMM4kR1gND0fthp21IfLkohDIYWBFK4YUB1MGxRVAgMF01HGMQp1DAAogiEqlKuDgiDOIUcBlQHSH7IKEWxh8EQg0AiDBA+AiEisQeBRBiEiDDAGIEQidjDYKhBIIUdBkBecHMIXMdoOcIgXqGHAdUB8qgQAhdzGAw5CCTCAPFJNhBSOMIo1jAYehBI4YcBME2ygRC7GMOAIIgH1QGmCW4OAYRB7EKvDggDzEKFEJjYwoAgOBFhgJgFHwhDumxFTGFAECwXehgAiwQfCEMQUxBIhME0MYTBEKuDFSev0toLz+27GdFIMhBiOsIopjAgCJaLIQikYYYByksyEGIRSxgQBNPFEgZAUQRCT0IPA0JgttiCgOoARREIPQg5DAiC2WILAqAsAqFjoYYBQTBfrGFAdYAyCIQOhRgGBMF8sQYBUAWB0JGQwoAQWCyFIKA6QFkEQgdCCQOCoJgUwgCogkBoWQhhQBAUk1IQUB2gCgKhRX2HAUFQTEpBANRBILSkzzAgCIpJNQioDlAVgdCCvsKAICgm1SAA6iIQGtZ1GBACxQwlBKgOUAeB0KAuw4AgKGYoQQA0gUBoCGEQDkIAqIZAaEBXYUAQzEYIAPURCDV1EQYEwXSEANAsAqEGwqAfBAHQDgKhorbDgCA4ESGwGEcYoS4CoYI2w4AgeAIhAHQryUC46JzTWruvMmHQLkIA6E+SgdCWtsJgyEFAAADhIBAKIgyaQwgAYSIQCiAM6iME0IeVJ5+ktRdu6LsZ0SAQ5mC+oB5CAIgLgTADYVANIQDEi0DoWGphQACEgXMQ0IRkA6HOoafMGcxGAISHMEBTkg2ERbhvQTEEADAcSQdC3/c0HospDAiAuFAdoElJB0IIYggDQiBOhAGaRiC0KNQwIAAATEMgtCC0ICAA0kN1gDYQCA0LIQwIgHQRBGgTgdCgvsKAAEgbIYCuEAgN6ToMCIH0EQToGoEQGYIgbYQA+kQgNKCL6oAgSBtBgBAQCDW1GQaEQNoIAYQm+EC4cN2p2nXwkb6bMVVbYUAQpINOHzEJPhBC1XQYEALxotNHKgiEkgiCYaLTxxAQCCU0GQYEQZjo+NEWM3u7pBdK2uful2TLflPS90l6VNK/SvpRd3+orzau6GvHsSEM0nHemtUzv4AW/aGkqyaW3SbpEnf/Rkn/LOn1XTcqL7gKYdPpq3XfoSN9N+METYUBQdAdOneExt3/xsw2Tyz7y9zDf5D0kk4bNSG4QAgNYRAmOnwk6MckvbfPBhAIczQRBgRBNXT4aIKdfJJWP3lTV7tbb2bbc4+3ufu2Iiua2RskPS7pXa20rCACYQbCoBt0/EjIAXffUnYlM9uq0WTzle7uzTerOAJhirphQBCciE4fmM7MrpL0Oknf6e69n4EbRSB0dbYyVUE9dPzAbGb2bknP1mhoabekN2l0VNFqSbeZmST9g7v/RF9tjCIQukAYFEOnD1Tj7i+fsvjGzhsyRzSB0GaVwBDRcnT8wPBEEwhS86FAVUDHD+AJUQWC1FwoDCkM6PTRtXWnrOy7CaggukCQwrgkdkhhQIePUBAEcYsyEKR6oRDrnAEdP0JGGMQv2kCQnujYywRDDGFAx4+YEATpKBwIZrZS0nZJ97n7C83sLI2uu7FZ0j2SXuruB9to5CJFq4U6YdBWEND5I1YEQXrKXP76NZJ25B5fJ+l2d79Y0u3Z494s6uxDCwMut4yYEQZpKhQIZna+pO+V9Lbc4qsl3ZT9fJOkFzXasgpmdfohhQFBgJitO2UlYZCwokNGb5X085LW5pZtcPc9kuTue8zs3IbbVsnkvEIIYUAAIHaEwDAsDAQzG9/y7U4ze3bZHZjZtZKulaSNmy4otE4TN8kJYfKYIEDMCIHhKVIhPEvS95vZCySdIul0M/tjSXvNbGNWHWyUtG/aytn1wLdJ0iWXXr7s0q7nrVmtB5bCukNanTAgBBArAgAL5xDc/fXufr67b5b0Mkl/5e6vkHSrpK3Zy7ZKuqW1VnaoahgwN4AYjecECANI9c5DuF7SzWb2Skm7JF3TTJP6UycMgBjQ8WOeUoHg7h+T9LHs5wclXdl8k/pBGCBVhACKivpM5aYQBkgJAYCqBh8IhAFSQAigCYMPhCoIA4SAEEDTBh0IZasDggB9ofNHFwYbCIQBQkTH3yxbdbJWbSh2QiwCDoQmzlZuCmGAJtHpI1TBBkKbylQHhAHmoXNHSgYXCIRBf+g8gbANLhCKCj0M6FwBNG1QgVC0Omg7DOjMAYRoMIHQZxgQAABiMJhAKKKJMKDzBxCrMvdU7lxTdyxr457Ik7iEMIDYBREIbY7ZdzFURBAASEEQgdA3wgAAEg+EtoeKCAMAKUk6EIqoWh0QBgBSM/hAqIIwAJCiZA87bWO4iCAAkLIkK4Q2jiwiDACkLslAaBphAGAIkguEpqsDwgDAUCQVCE3PGxAGAIYkmUAgDACgniQCoen7IxMGAIYoiUBoEmEAYKiiD4SmqwMAGKqoT0xj3gDAPLbqZK1Yf37fzYhGtBVCF/c4AIAhiTIQqoYBw0UAMFvwQ0abTl+t+w4dOeFxG5oYLlq3+sR8PXjkWO1tAkBXgg+EvDph0GZ1MBkEi5Z3gTACUFY0gdDmnEGd6qDPTn+eKu0iRIBhCyYQzluzWg8sHZn6XKgTyKGGQVWL/j0EBpC2YAKhTfOGi6pUB6kFQVHT/t2EBJCOQQRCk4YaBrMQEkA6kg+EpiaTCYLiZv2uCAogbMkHwjxFh4sIg2ZQTQBhG3QgFEEYtItqAghHUIEw70ijqtubZVF1QBD0i6AAuhdUIISCMAgXw05Ae5INhKrVAWEQH6oJDJ2ZnVXgZcfc/aF5L0gyELiIHSSuLYVBuT/7sjmvWSnpwnkbSS4Q6oQB1UHaGG5Cwna4++XzXmBmn1q0keACoc7EcpEw4CY4yKOKQCK+tYnXJPORuO4wEdUBpNH7gPcCYuPuX5EkM7vBzL5t3mvmCfKdX7ZzL/p6qgMUNQ4GwgGR+aSkXzSzfzGz3zSzLWVWjv7d3sQEMn/0mIdgQBPM7Ewze5+Zfd7MdphZkWGeUtz9Jnd/gaQrJP2zpF83sy8WXT+4OYQyyoQB1QHqyocCcw2o4AZJH3b3l5jZyZJObXFfT5P0tZI2S/pc0ZWC/dizqLPnonXoE0NKKMPMTpf0HZJulCR3f3TROQEV9zOuCH5F0t2Snunu31d0/SgrhLJh0HZ1sHJp/7JlR9ec0+o+EY5xKFA1YI6LJO2X9A4zu1TSnZJe4+4PN7yff5P0re5+oMrKQQdC09c2mlTn0920EJj3PAGRPoIhPL5ylY6tWd/V7tab2fbc423uvi37eZWkZ0h6tbvfYWY3SLpO0hub2LGZPUXSQ+7+v7PH3yXpRZL+XdLvuvujRbYTXb3bZ3Wwcmn/8a8u10VcGEoarAPuviX3tS333G5Ju939juzx+zQKiKbcLOk0STKzyyT9qaRdki6V9L+KbiToCkE6sUpo8pIURf9g2+jAqR6GgYoBY+7+gJnda2Zf4+5fkHSlSkz2FvAkd78/+/kVkt7u7r9tZisk3VV0I8EHwliVMKhaHXT9KX7e/giL+HF0EjKvlvSu7AijnZJ+tMFt569h9BxJr5ckdz9mNu/yRieKIhCaDoN51UFoQzpF2kNoxIOqYbjc/S5JpU4UK+GvzOxmSXskrZP0V5JkZhslFZo/kCIJhLKGds4BQ1DxIRjQsNdK+iFJGyV9u7s/li0/T9Ibim4kuUCocye00KqDqjgMNh4EA5rg7i7pPVOWL7zCaR6HQgxE/iinVIIvJRyVhDrM7M+beE1SFUKdoaKhdZL5fy/VQxjWrV5BpYCqvt3Mbp3zvEn6+kUbSSoQFuFT2HTjcCAYgGhdXeA1CyeXkwmEoU0kt4GqoX9UCajC3T/exHaSCYRFhjCZ3CTCoT+EAvqSxBgK1UG7mIjuHsOb6MNgKoRZ6OiKo2oAwmZmp0n6D3c/lj1eIekUd3+kyPqD+BjCp63mUTW0j/ctKrhdJ95451RJHy26cvQVAsNF/aJqaBfzCSjpFHdfGj9w9yUzK3xntuQ/gjCZ3B2qhnZQKaCEh83s+GW1zeyZkv6j6MpRVwhUB2GiagB681pJf2pm40thb9ToGkeFRB0IdfBJthuc9NYMho5QhLt/wsy+VtLXaHR28udzF7pbKNpAKFIdVC21Vywd6PK2e4NAMNRHKGAWM3uOu/+Vmf3gxFMXm5nc/QNFthNtILRlxdKBE74XQXgUx3BSPYQCZvhOje6B8H1TnnNJww6ELieTJ8ODgCiGqgFohru/Kfte6y5sUR6+EPpk8oqlA6UqjKHj6KRyOOoIs5jZ2Wb2P83sk2Z2p5ndYGZnF10/yXdW1eqg6U6cYCiHYCiOUMAM75G0X9KLJb0k+/m9RVeObsgo9OpgmnEoMJRUDENJaMyKVUN7H53l7r+ae/xmM3tR0ZWT+5hR58iitlEtlEPFAJT212b2MjNbkX29VNKHiq4cVSDUrQ5C6FwYRiovhP+3EDFshCl+XNKfaHQznEc1GkL6GTM7bGaHFq0c3ZDRPDH9gTCMVM7Kpf1DK/2B0tx9bZ31owmENquDPj+xcxJccYQCsJiZfb+k78gefszd/7zouvF8pF4gpupg0ngYiaGkxZhXOFHM73s0z8yul/QaSZ/Lvl6TLSskigohhbmDovKhQOUwG9UCMNULJF2Wu0HOTZI+Jem6Iisv/HhhZqeY2T+a2T+Z2d1m9svZ8rPM7DYz+2L2fV2Nf8RMbV6zSAr7yB8qh/liCnqgQ2fmfj6jzIpFetIjkp7j7pdKukzSVWb2LRolzu3ufrFGd+kplEBFrTtlZZTnHLSFcJiOUGDYCCf475I+ZWZ/mFUHd2bLClk4ZOTuLml8B56Tsi+XdLWkZ2fLb5L0MUmvK7rjWcqGwKI/hhQ7DI5QOhEnsgHH7598TNK3SPomjS5//Tp3f6DoNgp9tDCzlWZ2l6R9km5z9zskbXD3PZKUfT+3XPNPVKUiqPvJKPZP21QMJ0ox/IGisnmDV7n7Hne/1d1vKRMGUsFAcPej7n6ZpPMlXWFmlxTdgZlda2bbzWz7lx58cNnzbQ4NDaWDIBieMJT/80kMGyFzm5n9nJldkM3znmVmZxVdudRRRu7+kJl9TNJVkvaa2UZ332NmGzWqHqats03SNkm65NLLfby8bggMvTqYhqGkEY5AwoD9WPb9p3LLXNJFRVYucpTROWZ2ZvbzkyQ9V9LnJd0qaWv2sq2Sbimyw64mi4f6SVGiYpCG/f+PQfs6d39q/kvS1xdducjH7I0aXTDp05I+odEcwp9Lul7S88zsi5Kelz2ea+UKK9quuSiPixl6MBAKGKD/W3DZVEWOMvq0pMunLH9Q0pVFd9SUImGwqCMYWic55MtjDGn4iNtrDpeZnSdpk6QnmdnlGh1hJEmnSzq16HaiOFN5rIkwGCpCYRihgMH6bkk/otGBP2/JLT8s6ReKbiSqQFikSBgUqQ6OHditFevPb6JJQRnypDOhgJS5+02SbjKzF7v7+6tuJ5rB+CZOQCsaBvnvKRrakBkwILeb2VvGh/qb2W+bWeHLV0QRCH2djUwopIXhRAzAjRoNE700+zok6R1FVw4+EJoKgzLVwaJlqSAU0sMReIP3Ve7+JnffmX39sgqegyAFHgihvLkJBQCR+A8z+/bxAzN7lqT/KLpy1JPKbVcHk8+nONEsDW+ymQlmJOwnJL0zN29wUE+cQLxQsIEQ4lVMUw4FaViHphIKw/C4azDnZpjZSkmvcPdLzex0SXL3Q2W2EcaYzIQmw6DpIZGUh48khpCAWLn7UUnPzH4+VDYMpIArhFnaqAzKdvJUCmlItUrgjOVB+5SZ3SrpTyU9PF7o7h8osnJwgTCvOigbBm1+2iUU0pBqKGCwzpL0oKTn5Ja5pPgCoY8wqDMERCikgVBAKtz9R+usH0wgxFIZTCIUAITCzE6R9EpJT5d0yni5u//YzJVygpxU7kpTE8RDmGhOfbI59RPWMBh/JOk8jS5293GNLnZ3uOjKQQRCU9VBnx1X6qEgpX8EUkqhEMpJnejc09z9jZIezi54972SvqHoysm8a8p2Vm104IQCgJ49ln1/yMwukXSGpM1FVw46EIpewTSEMOhi22hfSlUCBmmbma2T9EaNbnP8OUm/XnTlYCaVywr5kyoTzQD64O5vy378uEpc1G6s9wph1ljnvE9qdcKgq0/wxw7sTrpaCDmQgaEys7PN7HfM7JNmdqeZvdXMzi66fu+BUEbdSeM+OmhCIT4MGyFi75G0T9KLJb1E0gFJ7y26cpBDRpN/kE10PH12zON9pziMxPAREJSz3P1Xc4/fbGYvKrpy8BVC7GGQF0o7sFgKVQKHng7SX5vZy8xsRfb1UkkfKrpyrxXCtDdsCn+I86RYLVAlAP0ys8MaXbPIJP2MRieoSdJKSUuS3lRkO0F/hEipOpgUaruqSnE+IfUPJ0iHu69199Oz7yvc/aTsa4W7n150O0EHQl2hd7qpHYmUYigAQxJUIOQ/kVXtXMadbEwdbUxtXSS1UKBKwJAEeZRRUSl1pCnNLTCnACxnZldJukGjcf23ufv1PTdpmd4CYXJCucgnsa4C4PG992rVhgs62dek1M9yjhH3S0Bd2f2Of0/S8yTtlvQJM7vV3T/X0PZPd/dDZnbWtOfd/UtFthNkhZAfdui6Cnh8773HvxMK1VElACe4QtK/uPtOSTKz90i6WqNrDTXhTyS9UNKdeuJoozFXwctYBBEIoYzTjsMgBIQCmsD9lYOxSVK+g9kt6Zub2ri7vzD7/tQ62wlqUlnqb1JyWhj0HRApzJGkMskcyocWBG29mW3PfV2be86mvN6bboCZPcvMTst+foWZvcXMLiy6fhAVwiwpdIh1pVApAH05esx18CtHu9rdAXffMuO53ZLyY9DnS7q/hTb8vqRLzexSST8v6UaNTlL7ziIrB1ch9GFeJdB3lSDFH4ypVAlADZ+QdLGZPdXMTpb0Mo3uV9C0x93dNZqfuMHdb5C0tujKBAI6kUIoMGyEqtz9cUmvkvQRSTsk3ezud7ewq8Nm9npJr5D0oezoppOKrtxLIHDRrfJirxKAoXP3v3D3r3b3r3L3X2tpNz8k6YikV7r7AxpNZv9m0ZV775lnnZ3c5TkHTbymC7GHQgpVAhAiMzvFzF4r6RclPSzp/0mSu+9y93cW3U7vgQDEhGEjBOomSVskfUbS90j67SobCfooI+nET+dNnygWyif/MmI/6ohzE4BWfL27f4MkmdmNkv6xykaiqhD67MBDCo/Yh44ANO6x8Q/ZBHYlQQbCvA6vqY45pA6+iphDIfa5BIaNEKBLzexQ9nVY0jeOfzazQ0U3EsyQUZlOos/rDAGx4fIV6XP3lU1sJ8gKYaytE8aqrhtaVUGVAKBJQQfCIo/vvTe4TrprhEI/GDZCinoNhGl/VFU6uDKhUDdAhh5AANIVdYWQN+SOmioBQBOCDYQqHfyidVIOjZhDAUAYggiEeZ8Sj9x/3/GvIrro9FMOlj7EWiUwj4DUBBEITRtqh02VAKCOoAJhUYdWtEpAfGKtEoCUBBUIY0P9hN8EqgQAVQUZCGPTKgKqBISEeQSkJOhAaFLTVUfIVUysVQLDRkC/ogwEqgQAaF4UgXB4196+mwAAyes8EMreT3lWGFAlzBfrsBGA/vR2+evxZNx43HjcgY3H5ic7/MO79mrthRs6bCH6wB3V2jHUS2A/dsz1wNKRvpsRjSiGjGahSpiPKqEbHGmEVEQVCKHNJYR8pBEAlBV0IIQWADGKsUrg8FOgH0EHwjSTITFr2IhP7wBQTpCBwNwAAHQvqECY9an+8K59E4+LVQkYYdgIQBFBBMK8DmsyDELD0BQkjjRCGoIIhGkmqwCqhHpirBIAdCvYQAAAdKvXQJg2TswnfYwxjwB0K6oKIfT5BACIWdCBsCgAOHGtHOYRAMwTTCBwtA4A9CuYQMjLf/I/uPMhHdz5UH+NAYCBCDIQ5glxHoHqBkAKoguEeThCaTHmEQDM0ksg5M/qLNpBzRo2YmIZAJoRVIWQ/4Q/b2goxGEjtINzEYDuBBUIQMy4nhFiF3QgTA4TcbQRALQnuECYNifw0L9/ecrrGDaqiollANMEEQjzDtucFgaTmFgGgPpW9d2AaeZ9+j+48yGtu+jM7hoDIFqPHj2m+w4d6bsZ0QiiQiiiSKXQJ05OAxC73gJh8nDCyZPKikwgM49QHfMIACZFUyEUxdnKAFDNwkAwswvM7K/NbIeZ3W1mr8mWn2Vmt5nZF7Pv69pvLgCgLUUqhMcl/ay7f52kb5H0U2b29ZKuk3S7u18s6fbscSeqnI+wasMFzTcEABKyMBDcfY+7fzL7+bCkHZI2Sbpa0k3Zy26S9KKmGxf6RDIApKTUHIKZbZZ0uaQ7JG1w9z3SKDQknVu3MYvOJyAgAKA9hQPBzNZIer+k17r7oRLrXWtm281s+4MHql2obNfhRyutBwAorlAgmNlJGoXBu9z9A9nivWa2MXt+o6Spx4C6+zZ33+LuW85ev750AwkDAOhGkaOMTNKNkna4+1tyT90qaWv281ZJtzTfPABAV4pcuuJZkn5Y0mfM7K5s2S9Iul7SzWb2Skm7JF1TdufTTo7iZDMA6MfCQHD3v5NkM56+stnmAAD6ktyZyvNwLgIAzBZkIEw78WzW5PJ4iIlLYANAPUEGwjzjcxG4exoANKv3QGjrstEMDwFAOb0HwiyclQwA3Qo2EAAA3Qo6EGI6S5khKgB1mNmrzewL2W0GfqOPNgR5T+U2rdpwAbe7BBAUM/suja4g/Y3ufsTMal8stIogKoRFdzm755HilcLqJ2+q25xBWLH+/L6bAOAJPynpenc/Iknu3sslG4IIBCAFR9ec03cTEK+vlvSfzOwOM/u4mX1TH40Y3JARALRkvZltzz3e5u7bxg/M7KOSzpuy3hs06ovXaXRXym/S6DpxF7m7t9ngScEHQpnhIgDIe/ToMe06+EhXuzvg7ltmPenuz531nJn9pKQPZAHwj2Z2TNJ6Sfubb+ZsDBkNEPMHQHD+TNJzJMnMvlrSyZKq3VGshuArBAzbsTXlb6oEROjtkt5uZp+V9KikrV0PF0kJBcLaCzf03QQAqMTdH5X0ir7bMcghI04iA4DlBhkITYspYJg/ADBLVIEQ06UsMCyhn4Nw8MixvpuACEQVCACA9hAICBZHGAHdSjoQGC8/Eb8PAPMkHQgAgOIIBACApIAC4fCuvZ3uL6ZDRQGgC8EEAmcat4v5AwCLBBMIVa29cPmNhfKf/tvuCKk02hHTEUahn4MAFBVVIFy49uS+mwAAyYoqEFBNbMNFMVUHQEoIBACApAgCYfOpDBPVEVt1EJsY5g+4jhGKCiIQVj9509znCYXhYLgI6E8QgVDFuovOXLZsUbA0LfQjjKgOAJQRdCA0dVTRrI4x9A59aKgOgH4FHQiojuqgfTHMHwBlJHNPZQCYdOSxY9q5/+G+mxGN6CqEM59yRt9NCF6M1QHDRUD/gg2EMh0/10ECgPp6DwQmdptFddCNWOYPOAcBZfQeCG3honYAUE40gVD3ENQYPzmXNYR/I4D2RBMISFOMw0VAqoIMhGlnIU8z7V4IQ0V10J1Y5g+AsoIMBJQTaxhQHQBhiTIQ5lUQZSd7q0wOhzShTBgAaEqvgTCtM2MYqLhYwyBmMQ0XccgpyoqqQph3slrXVzpFdVQHQJiCCoRpZxxfuPZk7qU8RazVQcxhEFN1AFQRVCBMKnL5ijKXrWiiEw1h/oAwANCGoANhnqHONRAG/aA6wBBEFwhFz1FIEWGAophQRhW9BcJkJzE5KTzZ8bd52euiw0B9DhfFGgYpoDrAUERXISAuVAdAPAiECMRaHaQQBlQHGJIgA2GoE8bTEAYoi/kDVBVkIISoj/kDwqBfVAcYmiACYV5ne+ZTzjg+oTyeaB5XEONzEMYT0kU67Vg62VjaCSAdQQRC02LvTGNuP9UBEK/gAqHMmcdNamJIaMX680/4GppUwgAYqlV9N2CedRedqYM7H+q7GZXDIh8Kxw7sLr1OTFIKg5irAyaUT/To48d0z76lvpsRjaADISWTHf20gIg1DACkofdAWLH+/OOd4+onb9KR+++TNJo4Prxr3/HXTU4oF3FszXqtWDowd5+LFK0OynbmqXT+VAZAOnqZQ2jiD2/yCKMmTHb+bYVBCo6tWZ9MGBxdc04SYcBwEeoKblJZOnFiuejF7EK4LPUQpBQEElUBkBdMIMzr0GedfzDNsrH6GZ3XrE/143ZQHTxhHAKpBUFKYUB1gCb0GgjTOpj8ENCi+YK2bptJGKQZAmMpBQHQpN4nlZtQ57DQopPL09ZNTYqdf16qQUB1gKYEMWQ0rXOdHBYqMlyUl+/cUu/o6ki5EshLNQyAJgURCGNlPuk3NVxU5ZN+7NXBUEJASm+uYBLVAZoU5JDR5PkIRXXRUccaBkPo/CelHARAG4IMhFmmDRcVrSpmnaQmFZ9LiC0MhhgC0nCCgOoATQs6ENZeuEGHd+1dtryto4tSMNQQkIYTBBJhgHYEHQh5fV0FdSzk6mDIISANKwiANvU+qTzuzMYd7ngIaFYVkF/e5NnJ8zr8EMNgSBPD8wwxDKgO0JbOA2H8Zi76h1y0MjjhUtMzOskqnWdIYUAIPCH1o4dmIQzQpiiGjLoaLqpzolpb6PxPNMQQALoSRSCMdT2ZXLY6mNZ5zzqyqex2ho4goDpA+4IOhPz5CJOqzh/MO/xUeqJKaCIM5i0/vr+sLYTAbIQBYYBu9D6pPE3Zzr7pcf6mwqDouoTBicbzA0OdJ5hEGKTPzK4xs7vN7JiZbcktf56Z3Wlmn8m+P6fNdvRaIRxdc45WLu0//ql93hj+rKOLil7uetprqgznlN0XpqOjX4wgGJTPSvpBSX8wsfyApO9z9/vN7BJJH5HU2th5cENG41BYteECPb733mXzBvOqh7Id9KLho6b2M1R0+tUQBMPj7jskycwml38q9/BuSaeY2Wp3P9JGO3oJhINHjmnd6tFo1WSVMM9kGCw61HS87XnKhAJBMBudf30EARZ4saRPtRUGUoFAMLO3S3qhpH3ufkm27CxJ75W0WdI9kl7q7gfrNmZy6GhcJUjVwiD/fV4wLAoFgmCETr89hEES1pvZ9tzjbe6+bfzAzD4q6bwp673B3W+Zt2Eze7qkX5f0/EZaOkORCuEPJf2upHfmll0n6XZ3v97Mrssev65qI+Z9ks+HwljRMCi6j/x28sGQchDQufePEGjf0ceP6aH9D3e1uwPuvmXWk+7+3CobNbPzJX1Q0n9x93+t2rgiFgaCu/+NmW2eWHy1pGdnP98k6WMqGQj5YaO8aRPMsyaRi4ZB/rkiQ0gpoMMPDwGAsszsTEkfkvR6d//7tvdXdQ5hg7vvkSR332NmxW9aMEORznrR4aBFOsEiQ0gho6OPCyGAIszsByT9jqRzJH3IzO5y9++W9CpJT5P0RjN7Y/by57v7vjba0fqkspldK+laSdp0/vzzC+Ydhrro8NKyHWWRACqzLUAiAFCNu39Qo2GhyeVvlvTmrtpRNRD2mtnGrDrYKGlmWmWTKtsk6dLLn+H552YNG+VNqwrqhsHkelWCgRCARAAgLVXPVL5V0tbs562S5s6QlzHuaIteCqKJjrnMNjh7dpgOHjk29QtISZHDTt+t0QTyejPbLelNkq6XdLOZvVLSLknXtNG4ycNB2wiD/LZmVQoEwLDQ0WOoihxl9PIZT13ZRAOmDRtN65zbDIPJba5c2k8IJIrOHpgtuEtXTCp7aGkTCIP40NED9QUbCLOGcIp21pMdxKLJa4SFDh7oXhCBMOtoo8lQqBoG+WUEQzvowIH4BREIRdQJg2nPtxEMRTvFNkOJjhlAVcEHQpnx/DKdYf61VTroOh0vnTaAEHUaCEeP+cznipykNm/dOqaFA502gKHpvEI4+JWjkqR1p6xc/lyJT+1tddgEAYCh6m3IaF4wSMvH+umoAaBdvc8hjINBWlw1AADaE9QxmAe/cvSEgGhiewCAYnqvEKZZNJyUf03RbcVo3r8fAJoWZCCMxdyZNyGUfz/BBAxD0IGAMHQdTAQQ0A8CAcFpK4AIGmA+AgGDUTdoCBSkLvhAeGDpyNTl561Z3cl+2tof4lMlUAiRfh19/Ji+/OAjfTcjGkEHwrxOetpzRTrtRR1/1fWa2Dehk54yIUJ4oG/BBkKVjrtqZ9+EJvZNYAxb0fAgONCWIANhsmO879ARbTqdzrDpwCNg4jQvOAgL1BFUIEwLgmk/Ew7NKBIwhEZcFlUZBAbmCSYQ5oXBpPFzTQfDrH0OOYDKViUESNhmBQZBASmQQFgUBrsOPqIL1526bL0mgmFe8Mx6zZADYpGiAUJwhIWggBRAIBQJg/z3JoKhSAiUWZ+AKI8J9Dgw0T0svQZC0TCYXDYtFPLrz+qg6wbBLG1sd+ghMy8wCIvwUGGkofcKoYp5oSAtPyqprSBo07Q2Dz0kxqqegwJgvmADYVp1MPn8vFAYKxsG+f0W2X6XCInZCAmgvmACYV7HvXP/w7ronNOWLS8yfFTGZAjNm7cIRdl/55AChJAAygkmEPLyHfPO/Q8f/z4tFNrY56znQg6GooZ+aC0hAcwWZCDMMi0Uig4dzbNoeGra61IIh7whD0dNhgQBgaEKOhDG1cEiVUOhaBDMWze1YMgbajVBQGCoVvS14/wfXb7jmdZJ37Nv6fjPs0KibOdeJwwmt9PUtmJx36Ejy75S9sDSkeNfQMqCrRDGHf84DO7Zt6TN5645/lzZSebJ1zWtqW3GWnEMpZqgekDKgg2EaYqEwiKLOu58BdLmJPYsk+2LNSDGUp+byAcE4YDYBRUIk51hfqgov2xeKMyqEop8gp8cjho/7iMYxlILCCndS3+Mw4FgQKx6D4RpnyAnO+aH9j+sM0t0ypOhUCUMpj3XZzCMTfu3xB4SqVURVA2IVe+BMMu4Ongo64zzoVBmPqFuGEx7XQjBkLfo3xhjYKRy/wuqBsQkmEAYd2rzOucqobBI0TCYtk5owTDLrN9DLEGRwhATVQNi0Nthp0WMq4MvP/jIsmVSscNR56myzuT6dbfRp/Ehs7EdOhv74a4cxopQ9VIhjP8Qxn/Qk9XBtMnkps3ryPPVR5ltxVIxzBLjBHbs1QNDSghJMENGkyargy8/+IjOOPvU48+VGTqatCgM8t+rBENVoQVK7AERUzgQDO04+tijWtp7b9/NiEYwgVCkOqgbCos67FmHuUrlgqGqvs+BWCS2gIgxHB5YOkIooDdBzyHk5w7atmiYqothrLwY5idimn+Iac6B+QX0JehAGMuXfF2GRAhCD4XYxBIKQB96D4Rpny4fKtgJljniiI61XTFUCQDm6z0QAEzHsBG6RiBEIJbqJpYqgWEjYLrgA4FDxjBkVAnoUvCBMFYnGGL5hI3uUCUAywUVCEUP7RzakUZSPKEWy7ARgOWCCoQqih6RBMSKYSN0JdhAGGIVkIpYqgSGjYATBRsIVXR9NnHXYhk2AhCnqAKhysQynWg/qBKaxbARuhBVIAAA2tNbIOQ/mVX5FL9ojiHVyiDVf1efqBKAkeAqhPxRQ+MhosN776m0LTrPfsUybARgJLhAqIJDTwGgvugCgUtZxFX5xFIlMGwERBgIbVl0R7Qu7phWVIh3U5sl9LuqASEws2vM7G4zO2ZmW3LLTzKzm8zsM2a2w8xe32Y7grmFZlFrNlywbNmZDXWQm89dc8K5DIRANbGFQCy315S453LCPivpByX9wcTyayStdvdvMLNTJX3OzN7t7ve00YggAuGic047Pgxy5jmnHZ8TWLPhAi3tvVdrN2xets743srztjlWZoiFEKgmthCQ4goCiTBImbvvkCQzW/aUpNPMbJWkJ0l6VNKhttrRWyBsOn114+O2szrzcccaw9g7IdC+2IJAIgwG7H2Srpa0R9Kpkn7a3b/U1s6CqBCmOePsU5eda1B3uChfiYQipgCQ4g0BiSBA69ab2fbc423uvm38wMw+Kum8Keu9wd1vmbHNKyQdlfRkSesk/a2ZfdTddzbV6LygAmFyDH+WRcNF8/RVLcTW8ecRAv0gDOo79tiRyucxVXDA3bfMetLdn1thm/9Z0ofd/TFJ+8zs7yVtkZR+IEwznkeYVh3k5YeLinS+bQRDzJ3+pJhDQIo7CCTCAMftkvQcM/tjjYaMvkXSW9vaWe+BcOG6U5cdq56fWJ6Urw6aOLqozjASARCe2INAIgyGyMx+QNLvSDpH0ofM7C53/25JvyfpHRodhWSS3uHun26rHb0HwjzjeYRF1UFdRauFlAJAIgRCRBgMk7t/UNIHpyxf0ujQ004EEwhlP6k3de7BZBukE4MhlRBIpfMfSykExggD9C2YQBibN7E8azK57PzBIrGHQGqd/1iKISARBAhHcIEwNm8eAU9ItfMfSzUExggDhCTYQBgbzyM0PZkcq9QDQEo/BMYIA4Qm+EBYJKRLTbSBAEgPQYBQRREIRU9EY+w/DkMLgDGCAKGLIhDyygwXTTvHIQRD6fjzhhoCEkGAeAQRCOOOe3zo6fhIo0UTy/OGi8ad7qzOt82gGGKHP2nIATBGECA2vQTCeWtW64GlI6WveNrkZPK0TntRSNDRT0fnfyKCALEKokKoYrI6aGL+gA5/MTr/2QgCxK73QBhXCYuGjfLVwbwwoFOvj05/MTp/pKi3QBgPG00zLRTG5p2VTBgsRmdfHp0/hqLXCmFyLmHeUUGLhoiGGgZ08M2i88eQ9T5kNCsU8lXCpCEMEdHRt4uOH1iu90CQFofCWJmqYLJDbfr+zUXQqfeHDh8or1YgmNlVkm6QtFLS29z9+kZapeXnJtQJg1nLEC86fKB5lQPBzFZqdDef50naLekTZnaru3+uyvbmnZtQdIiITj9+dPRAf+pUCFdI+hd33ylJZvYeSVdLqhQI0uJJ5llhQBCEjU4eiEOdQNgk6d7c492Svrlec2aHwrQwIAi6RccOpK1OINiUZb7sRWbXSro2e3jk6RvP+GyNfTZpvaQDfTciE1JbpLDaQ1umC6ktUjvteUrdDRx7eN9HHvn7t65vojEFhPT/UUmdQNgt6YLc4/Ml3T/5InffJmmbJJnZdnffUmOfjaEts4XUHtoyXUhtkcJrz5i7X9V3G2Kyosa6n5B0sZk91cxOlvQySbc20ywAQNcqVwju/riZvUrSRzQ67PTt7n53Yy0DAHSq1nkI7v4Xkv6ixCrb6uyvYbRltpDaQ1umC6ktUnjtQQXmvmweGAAwQHXmEAAACekkEMzsKjP7gpn9i5ld18U+F7TnHjP7jJndZWbbO973281sn5l9NrfsLDO7zcy+mH1f12NbfsnM7st+N3eZ2Qs6assFZvbXZrbDzO42s9dky/v63cxqT+e/HzM7xcz+0cz+KWvLL2fLO//dzGlLL+8bNKv1IaPsEhf/rNwlLiS9vOolLhpq0z2Strh758cNm9l3SFqS9E53vyRb9huSvuTu12eBuc7dX9dTW35J0pK7/1bb+59oy0ZJG939k2a2VtKdkl4k6UfUz+9mVnteqo5/P2Zmkk5z9yUzO0nS30l6jaQfVMe/mzltuUo9vG/QrC4qhOOXuHD3RyWNL3ExSO7+N5K+NLH4akk3ZT/fpFHH01dbeuHue9z9k9nPhyXt0Ohs+L5+N7Pa0zkfWcoenpR9uXr43cxpCxLQRSBMu8RFL39YOS7pL83szuxM6r5tcPc90qgjknRuz+15lZl9OhtS6mSIJs/MNku6XNIdCuB3M9EeqYffj5mtNLO7JO2TdJu79/a7mdEWqef3DerrIhAKXeKiY89y92dI+h5JP5UNnWDk9yV9laTLJO2R9Ntd7tzM1kh6v6TXuvuhLvddsD29/H7c/ai7X6bRFQGuMLNLuthvibb0+r5BM7oIhEKXuOiSu9+ffd8n6YMaDWv1aW82Zj0eu97XV0PcfW/2B39M0v9Rh7+bbEz6/ZLe5e4fyBb39ruZ1p4+fz/Z/h+S9DGNxux7fd/k29L37wXN6CIQgrrEhZmdlk0SysxOk/R8SX1fcO9WSVuzn7dKuqWvhow7mMwPqKPfTTZZeaOkHe7+ltxTvfxuZrWnj9+PmZ1jZmdmPz9J0nMlfV49/G5mtaWv9w2a1cmJadkhaG/VE5e4+LXWdzq7LRdpVBVIozO1/6TL9pjZuyU9W6OrQ+6V9CZJfybpZkkXStol6Rp3b32yd0Zbnq1R2e+S7pH04+Nx6pbb8u2S/lbSZyQdyxb/gkbj9n38bma15+Xq+PdjZt+o0aTxSo0+xN3s7r9iZmer49/NnLb8kXp436BZnKkMAJDEmcoAgAyBAACQRCAAADIEAgBAEoEAAMgQCAAASQQCACBDIAAAJEn/H9pI1/F12BowAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "psi = grid.cumint((ds.vo*ds.e3v).sum(dim='z_c'),'X', **bd) * 1e-6\n", "print(psi.coords)\n", "plt.contourf(\n", " domcfg.glamf,\n", " domcfg.gphif,\n", " psi.isel({'t':0}),\n", " levels=25,\n", " cmap='RdBu_r'\n", ")\n", "plt.colorbar(label='Psi barotropic [Sv]')\n", "plt.ylim(0,60)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By construction, $\\psi$ is 0 at the western boundary." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Kinetic Energy\n", "\n", "Finally, we plot the surface kinetic energy $1/2 (u^2 + v^2)$ by interpoloting both quantities the cell center point." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAJcCAYAAAAIM4iBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvjUlEQVR4nO3dfbRsdX3n+feHywUUZITgAwGMtk20MRPRxaBpTVrjQwPTHUzPShYmS+nEGUJ3GKUTJ2HSM2lnmekYWpNeyXKkrx06JGM0xoQO41xFYsdxTNTm6kIEkXAlqBcQBh940ATuPec7f9S+6ap99jm7at9Tdeoc3q+19jpV+/FX+9Q5v/rs32//KlWFJEkbOWqrCyBJWn5WFpKkXlYWkqReVhaSpF5WFpKkXlYWkqReVhaSpF5WFpKkXlYWHZLcleSVc9jvqUmuS3JPkkryzNbyY5NcneShJF9L8nOt5ZXk20keaaZ/31r+L5rtHmz2c2xPedZdf+wYh6eVJL+1wb5OTnJtU74vJ/mJsWXHJPlAc14ryct6ytV3Hs5O8pkk32l+nn0Er3Pdcq+zr1ck+WJz7D9L8j1jy5Lk15J8vZmuTJKN9rfMknwsyd+MvQdu71l/pveftpmqcmpNwF3AK+ew36cB/xz4AaCAZ7aW/yrw/wInAX8P+Bpw3tjyAv7uOvv+h8B9wPOa7T8GvG2Dsky9PnA88AjwQxvs773AHwAnAC8FHgSe1yw7Bri8mX8v8LKe87TueWj29WXgXwDHAm9snh8z5HVuVO6OfZ3SLP8x4Djg3wCfGlv+M8DtwOnAacAXgEu3+v18BO/XjwH//ZTrzvT+c9p+05YXYNkm4PeAVeCvm3+QvzCHYxy9TmVxN/DqsedvBd439nyjyuL3gX899vwVwNc2KMPU6wMXA3cCWWf58cBjwPe2zuOafxbAgSkqi3XPA/DqZnnGln+FsUp12tc5S7mbZZcAf9F63X8NPLd5/hfAJWPL3zBemfS85s8177fDU01xnl7aHPNbwFeBf7rJ79NZKouZ3n9O22/yMlRLVb2O0T+ff1xVJ1TVle11kjwjybc2mDa8lNElyUnAdzP6p3HY5xh9Uhv38Sbq/3HrMtbzOrZ9WpLvWueQs6x/MfC71fwX6PC9wEpV/WVP2XtNcR6eB9zcKsvNh5cneWmSb40t2+h19pa7+X2+tGtfVfVt4EutsvX9/jpV1fOb99sJwM8xSiifXW/9JM8APgT8FvAU4GzgpnXW/T82eK/e3FO0X03yQJI/77l8OOv7T9vM0VtdgO2oqr4CPHmTd3tC8/PBsXkPAk8ae/4PgE8BTwR+BfhgkrOr6lCzfXtbmu2/vs7xetdv/in9A0afkjcq+4Otee2yT6vvPGx4rKr6BJO/m41eZ2+5q6q9r/9vg/W7jnVCkmxQ0U5oKqZfAV5aVQ9tsOpPAn9aVe9tnn+d7t8zVfXPGV3+nNUvMrqU9hhwEfB/Ne+3L3WsO+v7T9uMyWJ5PNL8PHFs3onAw4efVNXHq+qxqvoW8CbgWYyu6R/evr0twMNJfnKskfJDfeu3yvV64BNV9VeHZyT50Nj+frJjX2vKPoO+8zDrsTZ6nUe6r76ynQg8MkNFcQbwfuDiVtrpcgajVDM3VfXpqnq4qh6tqmuAPwcuWGf1ad9P2qasLLpt+MfdXIZq9xZ6pPUPdLYDVn2TUePv88dmPx+4taech3vb3Nqx7X1V9fWqes/hSxxVdX7f+q1jvB64plXW88f29x7gL4Gjk5w5Q9m7X1D/ebgV+P5WL6Pv3+BYG73OWcs9sa8kxwPPbpVtlt/f30ryBOA/Av+2qj7UszqM2iiePeW+r9rgvTrL72j8/dY27ftJ29VWN5os48ToUs8lc9r3cYwaRgt4DnDc2LK3Af8Po94kz2X0T/NwL6DnMbouvYtR5P+3jK5r726Wn8eo19BZzfb/iY17Q/WuD/x94NvAk6Z4Xe9j1LPoeOAltHoVMeq5dByjBu5XN4/XazDf6Dwc7g31pmafl7Fxb6gNX2dfuVv7ekqz/L9ryv9rTPaGuhS4jVFPqO9m9A/00rHlHwPess6+3wu8Z4b30TMYfWr/cUaXk78LOHsT36dPZtTD6bhm/z/ZvBeeM/T95LS9py0vwDJOwIWMGrm/Bbx5k/dd7Wls2bHA1cBDjLoh/tzYsh9uKodvA/cz+hR6ZmvfP9ds9xDwH4Bje8qy4frAvwN+b8rXdXJTpm835+4nWsvv6njtz1xnX+ueh2b5C4DPMOqJ9FngBWPLfpDRpZ+pXucU5X4E+MGx568Evtgc+2Pjr4HRp+4rgW8005VM9tr6EvCqDd4X32GyR9QPdq3beq2fbl7XVxldvtqs9+lTgBsZVUjfYvQB6lVjy5/RlPEZQ99/TttrSvNLljRHSU4H/rCqfmCryyINYWUhSeplA7e0pFq92IY2SkubwmQhSeq1LW7KOybH1nEcv9XFkLQNPMw3H6iqpxzJPv7hy4+vr39jZbOKtKHP3Pzo9VV13kIOdgS2RWVxHMfzorxiq4shaRv40/rAl490H1//xgr/+fpnbEZxeu069Y5TFnKgI2SbhSSp17ZIFpK0SAWssrrVxVgqJgtJUi+ThSStUayUyWKcyUKS1MtkIUktozYL70EbZ7KQJPUyWUhSB3tDTTJZSJJ6mSwkqaUoVhw3b4LJQpLUy2QhSR3sDTXJZCFJ6mVlIUnq5WUoSWopYMXLUBNMFpKkXiYLSepgA/ckk4UkqZfJQpJaCrwpr8VkIUnqZbKQpA4OIzjJZCFJ6mWykKSWorzPosVkIUnqZbKQpLaCFYPFBJOFJKmXyUKSWgp7Q7WZLCRJvUwWkrRGWCFbXYilYrKQJPWyspAk9fIylCS1FLBq19kJJgtJUi+ThSR1sIF7kslCktTLZCFJLYXJos1kIUnqZbKQpA6rZbIYZ7KQJPUyWUhSi20Wa5ksJEm9TBaS1FKEFT9LT/BsSJJ6mSwkqYO9oSaZLCRJvUwWktRib6i1TBaSpF5WFpKkXl6GkqQ1wkr5WXqcZ0OS1MtkIUktBaz6WXqCZ0OS1MtkIUkd7Do7yWQhSeplspCklip7Q7V5NiRJveZWWSQ5Lsl/TvK5JLcm+d+a+W9JcneSm5rpgnmVQZKGWiULmbaLeV6GehT44ap6JMlu4BNJPtQs+42qevscjy1J2kRzqyyqqoBHmqe7m6nmdTxJ2iyjgQS9Sj9urmcjya4kNwH3AzdU1aebRZcluTnJ1UlOWmfbS5LsS7LvII/Os5iSpB5zrSyqaqWqzgZOB85N8n3Au4BnA2cD9wLvWGfbPVV1TlWds5tj51lMSWoZ9YZaxNRbkuS8JLcn2Z/kio7lz03yySSPJnnz2PznjLUN35TkoSSXN8tmbjteSNfZqvpWko8B5423VSR5N/DBRZRBkrabJLuAdwKvAg4ANya5rqq+MLbaN4A3Aq8Z37aqbmf0ofzwfu4Grh1bZaa243n2hnpKkic3j58AvBL4YpJTx1b7UeCWeZVBkoY4PDbUIqYe5wL7q+rOqnoMeB9w4URZq+6vqhuBgxvs5xXAl6rqy0PPyTyTxanANU2NdhTw/qr6YJLfS3I2o9/HXcDPzLEMkrTsTkmyb+z5nqra0zw+Dfjq2LIDwIsGHOMi4L2teZcleT2wD/j5qvrmRjuYZ2+om4EXdMx/3byOKUnb0ANVdc46y7puxJipV2mSY4AfAf7nsdnvAt7a7OutjNqOf3qj/TjchyR1WKmluGHuAHDG2PPTgXtm3Mf5wGer6r7DM8YfT9t2bEdiSVpeNwJnJnlWkxAuAq6bcR+vpXUJakjbsclCklqKLMVNeVV1KMllwPXALuDqqro1yaXN8quSPJ1Ru8OJwGrTPfasqnooyRMZ9aRqtw1fOWvbsZWFJC2xqtoL7G3Nu2rs8dcYXZ7q2vY7wHd1zJ+57djKQpI6rDpE+QTPhiSpl8lCklocSHAtz4YkqZfJQpJaiizLfRZLw2QhSeplspCkDlMM8ve44tmQJPUyWUhSSxVTfTHR44lnQ5LUy2QhSWuE1c7RwR+/rCwkzV8GXMSo1c0vhwbzMpQkqZfJQpJaChu42zwbkqReJgtJ6uBAgpM8G5KkXiYLSWopwqoDCU4wWUiSepksJKmDbRaTPBuSpF4mC0lqKWDV+ywmeDYkSb1MFpK0RlhxIMEJJgtJUi+ThSS12GaxlmdDktTLZCFJHWyzmGSykCT1MllIUktVbLNo8WxIknpZWUiSenkZSpI6+LWqkzwbkqReJgtJailg1a6zE0wWkqReJgtJWiO2WbR4NiRJvUwWktQyGkjQNotxJgtJUi+ThSR1WPGz9ATPhiSpl8lCklqK2GbRYrKQJPUyWUhSh1U/S0/wbEiSepksJKmlClZss5hgspAk9bKykCT18jKUJHWw6+wkk4UkqZfJQpJaRjfl+Vl6nGdDktTLZCFJHVb8WtUJJgtJUq+5JYskxwEfB45tjvOBqvpXSU4G/gB4JnAX8ONV9c15lUOSZuWXH601z2TxKPDDVfV84GzgvCQvBq4APlpVZwIfbZ5LkpbY3JJFVRXwSPN0dzMVcCHwsmb+NcDHgF+cVzkkaXb2hmqb69lIsivJTcD9wA1V9WngaVV1L0Dz86nrbHtJkn1J9h3k0XkWU5LUY66VRVWtVNXZwOnAuUm+b4Zt91TVOVV1zm6OnVsZJanLKlnI1CfJeUluT7I/yZrL9kmem+STSR5N8ubWsruSfD7JTUn2jc0/OckNSe5ofp7UV46F5Kyq+hajy03nAfclORWg+Xn/IsogSdtNkl3AO4HzgbOA1yY5q7XaN4A3Am9fZzcvr6qzq+qcsXkztx3PrbJI8pQkT24ePwF4JfBF4Drg4ma1i4E/mVcZJGmIw0OUL2LqcS6wv6rurKrHgPcxavcdK2vdX1U3AgdneIkXMmozpvn5mr4N5nlT3qnANU3NeBTw/qr6YJJPAu9P8gbgK8CPzbEMkrTsThm/RATsqao9zePTgK+OLTsAvGiGfRfwkSQF/Lux/U60HSfpbDseN8/eUDcDL+iY/3XgFfM6riRthgX2hnqgdYloXFf0qBn2/ZKquqepDG5I8sWq+vjsRfQObklaZgeAM8aenw7cM+3GVXVP8/N+4FpGl7VgQNuxlYUkLa8bgTOTPCvJMcBFjNp9eyU5PsmTDj8GXg3c0iyeue3YgQQlTS+Pj8+XoyHKt364j6o6lOQy4HpgF3B1Vd2a5NJm+VVJng7sA04EVpNczqjn1CnAtUlg9L/+96vqw82u38aMbcdWFpK0xKpqL7C3Ne+qscdfY3R5qu0h4Pnr7HPmtmMrC0nqMM0Nc48nj49MKUk6IiYLSWpxiPK1TBaSpF4mC0nq4BDlkzwbkqReJgtJaqvluM9imZgsJEm9TBaS1FJ4n0WbyUKS1MtkIUkdbLOYZGUhae5y1Oz/eGtlDgXRYFYWktTiHdxr2WYhSeplZSFJ6uVlKEnq4GWoSSYLSVIvk4UktSzL16ouE5OFJKmXyUKSOjjcxySThSSpl8lCktrK3lBtJgtJUi+ThSS1ONzHWiYLSVIvk4UkdTBZTDJZSJJ6mSwkqcU7uNcyWUiSepksJKlDmSwmmCwkSb2sLCRJvbwMJUkdHEhwkslCktTLZCFJLeVAgmuYLCRJvUwWktTBrrOTTBaSpF4mC0law+E+2kwWkqReJgtJ6mCbxSSThSSpl8lCklr8WtW1TBaSpF4mC0lqq9Fd3PovTBaSpF4mC0nq4Kizk0wWkqReVhaSpF5ehpKklsKb8tpMFpKkXiYLSVrDgQTbdm5lkQWGplpd3LEkaQvs3MpCko6AN+VNss1CktRrbpVFkjOS/FmS25LcmuRNzfy3JLk7yU3NdMG8yiBJQ1VlIdN2Mc/LUIeAn6+qzyZ5EvCZJDc0y36jqt4+x2NLkjbR3CqLqroXuLd5/HCS24DT5nU8SdosVd5n0baQNoskzwReAHy6mXVZkpuTXJ3kpHW2uSTJviT7DvLoIoopSVrH3CuLJCcAfwRcXlUPAe8Cng2czSh5vKNru6raU1XnVNU5uzl23sWUpAmrlYVM28VcK4skuxlVFO+pqj8GqKr7qmqlqlaBdwPnzrMMkrSdJTkvye1J9ie5omP5c5N8MsmjSd48Nr+zk1GzbOaORnNrs0gS4LeB26rq18fmn9q0ZwD8KHDLvMogSUMtw30WSXYB7wReBRwAbkxyXVV9YWy1bwBvBF7T2ryzk9HYtjN1NJpnb6iXAK8DPp/kpmbeLwGvTXI2o7G67gJ+Zo5lkKTt7Fxgf1XdCZDkfcCFwN9WFlV1P3B/kv92fMMNOhmNVzRTm2dvqE9A57eH7J3XMSVpsyywN9QpSfaNPd9TVXuax6cBXx1bdgB40awH6OhkBKOORq8H9jFKIN/caB/ewS1JW+uBw515mmnP2LKuGmumC2QdnYxgyo5G46wsJGl5HQDOGHt+OnDPtBt3dTKCYR2NduxAgtm1a+ZtamVl4MEWVOc6uu2RG/K78rw/7hRLMxTHjcCZSZ4F3A1cBPzENBuu18moWTZzR6MdW1lI0nZXVYeSXAZcD+wCrq6qW5Nc2iy/KsnTGbU7nAisJrkcOAv4fjo6GVXVXuDKWTsaWVlIUocl6DkLQPPPfW9r3lVjj7/G6PJU23qdjKiq181aDtssJEm9TBaS1OZAgmuYLCRJvUwWktRlWRotloTJQpLUy2QhSR1ss5hkspAk9TJZSFKHZRiifJmYLCRJvUwWktRS2GbRtnMri0UO/jbkWIsafHDRx1r2QfeWvXxLLkf5D/TxaudWFpI0VAEmiwm2WUiSellZSJJ6eRlKkjrYdXaSyUKS1MtkIUldTBYTTBaSpF4mC0laI96U12KykCT1MllIUhfbLCaYLCRJvUwWktRWDiTYZrKQJPXaucliwEirOWrYiKS1Mmiz5TZ0dNYhI9zu1JFgPRfbm20WE0wWkqReOzdZSNIRsc1inMlCktTLZCFJXWyzmGCykCT1srKQJPXyMpQkdfEy1ASThSSpl8lCktoKcLiPCSYLSVIvk4UkdSjbLCaYLCRJvUwW44YM/AbA7CMJZteumbepQ8MGmctRs197HTo44iKPpcbQ9+0CBy2s1W34MX0bFnmeTBaSpF4mC0nqYm+oCSYLSVIvk4UkdYhtFhNMFpKkXiYLSWor7A3VYrKQJPUyWUjSGrE3VIvJQpLUy8pCktTLy1CS1MUG7gkmC0lSL5OFJHUxWUzYsZXFkNFPFzl657BRZw/OvM3oYENe1+BhZxd3rGW3wFFdpXnbsZWFJB0Rk8UE2ywkSb3mVlkkOSPJnyW5LcmtSd7UzD85yQ1J7mh+njSvMkjSIMXoprxFTNvEPJPFIeDnq+rvAS8GfjbJWcAVwEer6kzgo81zSdISm1ubRVXdC9zbPH44yW3AacCFwMua1a4BPgb84rzKIUlDOET5pIW0WSR5JvAC4NPA05qK5HCF8tR1trkkyb4k+w7y6CKKKUlax9x7QyU5Afgj4PKqeiiZ7hpdVe0B9gCcmJOt4yUtlv91JvQmiyT/OsmTx56flORXptl5kt2MKor3VNUfN7PvS3Jqs/xU4P6ZSy1JjxNJzktye5L9Sda08SZ5bpJPJnk0yZun2XZIR6NpLkOdX1XfOvykqr4JXNC3UUYR4reB26rq18cWXQdc3Dy+GPiTKcogSY87SXYB7wTOB84CXtt0FBr3DeCNwNtn2HbmjkbTVBa7khw7VoAnAMdusP5hLwFeB/xwkpua6QLgbcCrktwBvKp5Lkla61xgf1XdWVWPAe9j1Enob1XV/VV1I9Ae4mGjbS9k1MGI5udr+goyTZvF/wl8NMl/YHQV76fHDrKuqvoEsF4DxSumOK4kbZkF9oY6Jcm+sed7mjZbGPUg/erYsgPAi6bc70bbTnQ0StLZ0Whcb2VRVVcmuRl4JaN//m+tquunLKwkaWMPVNU56yzr+sA9bTV2JNuuMVVvqKr6MPDhztIkn6yqHxhagLk5evaOXjlqWE/i1ccem/1YDjInLbfluLv6AHDG2PPTgXs2Ydv7kpzapIqpOhptxn0Wx23CPiRJa90InJnkWUmOAS5i1EnoSLeduaPRZtxnYW9kSZqDqjqU5DLgemAXcHVV3Zrk0mb5VUmeDuwDTgRWk1wOnNXc17Zm22bXbwPen+QNwFeAH+sri0OUS1JbsTQfg6tqL7C3Ne+qscdfY3SJaaptm/lfZ8aORtPclHdZzw0bS3FhT5I0P9O0WTwduDHJ+5u7AduVw+vmUC5J2lq1oGmb6K0squp/Ac5kdDf2PwXuaIYAeXaz/Ja5llCStOWm6g1VVQV8rZkOAScBH0hy5RzLJklbJrWYabvobeBO8kZGXaseAP498D9V1cEkRwF3AL8w3yJKkrbaNL2hTgH+SVV9eXxmVa0m+UfzKZYkbbFt9Kl/EaYZ7uOXN1h22+YWR5K0jLzPQpK6mCwmLORrVSVJ25vJQpJatltPpUXYuZVFzf6brtUlHwk2w4Jgjpr9JvuhfyeDjjXwdS09RxbWDrJzKwtJOhLLMUT50tihH+kkSZvJZCFJXWyzmGCykCT1srKQJPXyMpQkdbDr7CSThSSpl8lCkrqYLCaYLCRJvUwWktTmcB9rmCwkSb1MFpLUxWQxwWQhSeq1Y5NFdu2afaNjdg871qOPDthogfX0Ao9Vq7N/HBs0Uu2A42gTDH4vrWxqMRbCt9gEk4UkqdeOTRaSdCTsDTXJZCFJ6mVlIUnqZWUhSeplm4UkdbHNYoLJQpLUy8pCktTLy1CS1OZAgmuYLCRJvUwWktTFZDHBZCFJ6mWyGJOBg6QNGjyvVgcda5ABA/UNHTBuyACOdWjAuRh6/hY5gKO2N5PFBP9yJEm9TBaS1BLsDdVmspAk9TJZSFIXk8UEk4UkqZfJQpLavIN7DZOFJKmXyUKSupgsJpgsJEm9TBaS1MVkMcFkIUnqZWUhSerlZShJ6mDX2Uk7t7IYMvrpwJFMM2RU1wHl4+Ch2bcBGDIq7pDyAewaEFZXBmzjX7K0UDu3spCkI+HnkQm2WUiSes2tskhydZL7k9wyNu8tSe5OclMzXTCv40vSYLXAaZuYZ7L4HeC8jvm/UVVnN9PeOR5fkrRJ5tZmUVUfT/LMee1fkubJPhSTtqLN4rIkNzeXqU5ab6UklyTZl2TfQR5dZPkkSS2LrizeBTwbOBu4F3jHeitW1Z6qOqeqztnNsQsqniQ1bLOYsNDKoqruq6qVGt3Q8G7g3EUeX5K2myTnJbk9yf4kV3QsT5LfbJbfnOSFzfznjHUmuinJQ0kub5bN3NloofdZJDm1qu5tnv4ocMtG60vSVlmGNosku4B3Aq8CDgA3Jrmuqr4wttr5wJnN9CJGV3BeVFW3M7qKc3g/dwPXjm33G1X19mnLMrfKIsl7gZcBpyQ5APwr4GVJzmYUvu4CfmZex5ekHeBcYH9V3QmQ5H3AhcB4ZXEh8LtVVcCnkjy59cEc4BXAl6rqy0MLMs/eUK/tmP3b8zqeJG2qxSWLU5LsG3u+p6r2NI9PA746tuwAo/Qwrmud0xi1Cx92EfDe1naXJXk9sA/4+ar65kaF9A5uSdpaDxzuzNNMe8aWdQ08167GNlwnyTHAjwB/OLZ86s5Ghzk21JhkWN1ZQwbqG7DNULWysrBjDTJwAMcdach7cJHnb+ixttvveHl6Kh0Azhh7fjpwz4zrnA98tqruOzxj/HGSdwMf7CuIyUKSlteNwJlJntUkhIuA61rrXAe8vukV9WLgwVZ7xWtpXYJKcurY06k6G5ksJGlJVdWhJJcB1wO7gKur6tYklzbLrwL2AhcA+4HvAD91ePskT2TUk6rdmejKWTsbWVlIUkvobgjYCs0Yentb864ae1zAz66z7XeA7+qY/7pZy+FlKElSL5OFJHVZjgbupWGykCT1MllIUodlGO5jmZgsJEm9TBaS1MVkMcFkIUnqZbKQpC4miwkmC0lSL5OFJLWVvaHarCw2wVG7B5zGo2YfTCADtgHIrl2zbzR0lNCVxYwuOvRcDBkhWP+F5+/xy8pCkrpYL06wzUKS1MtkIUkdbLOYZLKQJPWyspAk9fIylCR18TLUBJOFJKmXyUKSOtjAPclkIUnqZbKQpLbCNosWk4UkqZfJQpK6mCwmmCwkSb12bLLI0Yt7aUNG4swiu1oMGUE2i/scMej8DRx1dtnPxdIbOhrxNhPsDdXmX4EkqdeOTRaSdERMFhNMFpKkXiYLSeqQMlqMM1lIknqZLCSpzTu41zBZSJJ6WVlIknp5GUqSOnhT3iSThSSpl8lCkrqYLCaYLCRJvXZssqghA54dPDToWIMGtVvggGxDBuqDlU0vx3oGnb+hg/sNuRC9yMHzHicD9W0HtllMMllIknrt2GQhSUfEZDHBZCFJ6mWykKS2ss2izWQhSeplspCkLiaLCSYLSVIvk4UktQTbLNpMFpKkXiYLSeri16pOMFlIknpZWUiSenkZSpI62MA9acdWFjl69pdWK0s+EuzQkVaHjGS6yGMNOczK4kbFlbSDKwtJGqzwprwW2ywkSb3mVlkkuTrJ/UluGZt3cpIbktzR/DxpXseXpCOR1cVM28U8k8XvAOe15l0BfLSqzgQ+2jyXJC25uVUWVfVx4But2RcC1zSPrwFeM6/jS9IRqQVN28Si2yyeVlX3AjQ/n7reikkuSbIvyb6DPLqwAkqS1lra3lBVtQfYA3BiTt5G9a+kncD7LCYtOlncl+RUgObn/Qs+viRpgEVXFtcBFzePLwb+ZMHHl6R+xWggwUVMPZKcl+T2JPuTrOkUlJHfbJbfnOSFY8vuSvL5JDcl2Tc2f+aeqfPsOvte4JPAc5IcSPIG4G3Aq5LcAbyqeS5J6pBkF/BO4HzgLOC1Sc5qrXY+cGYzXQK8q7X85VV1dlWdMzZv5p6pc2uzqKrXrrPoFfM6piRtliVpszgX2F9VdwIkeR+jXqVfGFvnQuB3q6qATyV5cpJTD3cmWseFwMuax9cAHwN+caOCeAe3JG2tUw73/GymS8aWnQZ8dez5gWYeU65TwEeSfKa136l7ph62tL2hjtiALy6pgweHHWvXrtm3WT008yY5KrMfB6gFjrk3ZIDEDDl/Nex3NeRYizx/WiKLSxYPtC4Rjev6o2+XbKN1XlJV9yR5KnBDki8298DNzGQhScvrAHDG2PPTgXumXaeqDv+8H7iW0WUtGNAz1cpCkpbXjcCZSZ6V5BjgIka9SsddB7y+6RX1YuDBqro3yfFJngSQ5Hjg1cAtY9vM1DN1516GkqSBwnI0cFfVoSSXAdcDu4Crq+rWJJc2y68C9gIXAPuB7wA/1Wz+NODaJDD6X//7VfXhZtnbgPc3vVS/AvxYX1msLCRpiVXVXkYVwvi8q8YeF/CzHdvdCTx/nX1+nRl7plpZSFLblDfMPZ7YZiFJ6mWykKQOy9BmsUxMFpKkXiYLSepisphgspAk9TJZSFIH2ywmmSwkSb1MFpLUVsCAQTF3sh1bWdTK7EOF1mMDRzLdPftpHFK+0d3+sxtyrAw7FNTqgG2GjaYraXF2bGUhSUfEYDHBNgtJUi+ThSR1sDfUJJOFJKmXlYUkqZeXoSSpi0OUTzBZSJJ6mSwkqYMN3JNMFpKkXiYLSWorvCmvxWQhSeplspCklgCxN9QEk4UkqdfOTRZDRpAdMmIqUI89tpBj1aFh5Rti2Ki4A4+1wKGgF3ksbXOL+3PbFkwWkqReOzdZSNIRsM1ikslCktTLZCFJbd5nsYbJQpLUy2QhSWuUo862mCwkSb1MFpLUwVFnJ5ksJEm9rCwkSb28DCVJXWzgnmCykCT12rHJolZnHwXsqBNOGHasv/7r2Y/1hCfMfpwhAxYC7No1+zZLP5DggNc0OtjCjrX0gxYOHDjzcaEgnp4JJgtJUq8dmywk6YjYZjHBZCFJ6mWykKQuBosJJgtJUi+ThSR18MuPJpksJEm9TBaS1MVkMcFkIUnqZbKQpLYCvIN7gslCktTLZCFJLaHsDdVispAk9dqxyeLQuc+deZuDJww7Hcff8Y2Zt1l94jEzb3PUN7898zYAHD1g1NS/eXTYsYY4eHBhh6oFHouDh2bfZsCn2RpynIFq4GjEg7ZzVNylsmMrC0k6Il6GmuBlKElSry1JFknuAh4GVoBDVXXOVpRDktZlspiwlZehXl5VD2zh8SVJU7LNQpLavClvja1qsyjgI0k+k+SSrhWSXJJkX5J9B1lgzxxJ0hpblSxeUlX3JHkqcEOSL1bVx8dXqKo9wB6AE3OyFw8lLZQ35U3akmRRVfc0P+8HrgXO3YpySJKms/DKIsnxSZ50+DHwauCWRZdDkjZUtZipR5LzktyeZH+SKzqWJ8lvNstvTvLCZv4ZSf4syW1Jbk3yprFt3pLk7iQ3NdMFfeXYistQTwOuTXL4+L9fVR/egnJI0lJLsgt4J/Aq4ABwY5LrquoLY6udD5zZTC8C3tX8PAT8fFV9tvmA/pkkN4xt+xtV9fZpy7LwyqKq7gSev+jjStL0pvvUvwDnAvub/5skeR9wITBeWVwI/G5VFfCpJE9OcmpV3QvcC1BVDye5DTitte3UvINbkrbWKYd7fjbTeA/R04Cvjj0/0MxjlnWSPBN4AfDpsdmXNZetrk5yUl8hd+x9Fn/1P8z+qeC4W4adjifcPfuggCsnHDvzNnlscQPGZcjgg0CNLi/OdqwFDbg3OtaAgQRXBna43z3gdR014PPbowO7lg8574dm//0C5KjZt1t97LFBx9oUxSKTxQMbjGLRdeLaBdtwnSQnAH8EXF5VDzWz3wW8tVnvrcA7gJ/eqJAmC0laXgeAM8aenw7cM+06SXYzqijeU1V/fHiFqrqvqlaqahV4N1P0SLWykKQuqwuaNnYjcGaSZyU5BrgIuK61znXA65teUS8GHqyqezPqRfTbwG1V9evjGyQ5dezpjzJFj9QdexlKkra7qjqU5DLgemAXcHVV3Zrk0mb5VcBe4AJgP/Ad4KeazV8CvA74fJKbmnm/VFV7gSuTnM3oMtRdwM/0lcXKQpKWWPPPfW9r3lVjjwv42Y7tPkF3ewZV9bpZy2FlIUkdHO5jkm0WkqReJgtJ6mKymGCykCT1MllIUlsBqyaLcSYLSVIvk4UkrbE0AwkuDZOFJKmXyUKSupgsJuzYymL/y39n5m2e+1f/bNjBvnDnzJvs2jV7qFt9bMCIqQwc8XNlZdCxyM4Lq9k1bATeQQb8rgY3xA451oBRhQE4esC/mq0cdVZr7NjKQpKOiMliws77GChJ2nQmC0lq8z6LNUwWkqReJgtJWqOgBn6d7g5lspAk9bKykCT18jKUJHWx6+wEk4UkqZfJQpLa7Dq7hslCktTLZCFJXWyzmGCykCT12rHJ4s//ZvYbalZ3D/wkMWAE2ezePexYC5Kho8cOuJGphlwbXuANUzV4AN4BI7QOOdbQ39WAEWQXOgLvVjNZTDBZSJJ67dhkIUnD+bWqbSYLSVIvk4UktRWw6kCC40wWkqReJgtJ6mKbxQSThSSpl8lCkrqYLCaYLCRJvawsJEm9vAwlSWuUQ5S3mCwkSb12bLJ43f/9z2be5okPDhj4DVj967+ZfaMh2wy1wEH3Bg9qN+thBg5ol92zv+Vz9MA/kwFlHFI+jj129m0Ajpl9MMt6wjGDDpVvPTL7Rg8/POhYm6KgFvl3sw2YLCRJvXZsspCkI2KbxQSThSSpl8lCkrp4U94Ek4UkqZfJQpLaqhyivMVkIUnqZbKQpC62WUwwWUiSepksJKlD2WYxwWQhSeplspCkNco2ixaThSSp145NFmf+j5/a6iI8Pi1opM46NOw4dejgJpdEenzYsZWFJA1WOJBgi5ehJEm9tqSySHJektuT7E9yxVaUQZI2VKuLmbaJhVcWSXYB7wTOB84CXpvkrEWXQ5I0va1oszgX2F9VdwIkeR9wIfCFLSiLJK1RQNlmMWErLkOdBnx17PmBZt6EJJck2Zdk30EeXVjhJElrbUWySMe8NVV4Ve0B9gCcmJOt4iUtTtW2ak9YhK1IFgeAM8aenw7cswXlkCRNaSuSxY3AmUmeBdwNXAT8xBaUQ5LWZZvFpIUni6o6BFwGXA/cBry/qm5ddDkkaTvou9UgI7/ZLL85yQv7tk1ycpIbktzR/Dyprxxbcp9FVe2tqu+tqmdX1f++FWWQpA0twX0WU95qcD5wZjNdArxrim2vAD5aVWcCH22eb8g7uCVpef3trQZV9Rhw+FaDcRcCv1sjnwKenOTUnm0vBK5pHl8DvKavINtibKiH+eYDf1of+PI6i08BHlhkeZa0DLAc5ViGMsBylGMZygDLUY5FluF7jnQHD/PN6/+0PnDKZhRmCscl2Tf2fE/TGxS6bzV4UWv79W5H2Gjbp1XVvQBVdW+Sp/YVcltUFlX1lPWWJdlXVecssjzLWIZlKccylGFZyrEMZViWcixDGWZRVedtdRka09xqsN46U92mMC0vQ0nS8prmVoP11tlo2/uaS1U0P+/vK4iVhSQtr7+91SDJMYxuNbiutc51wOubXlEvBh5sLjFttO11wMXN44uBP+kryLa4DNVjT/8qc7cMZYDlKMcylAGWoxzLUAZYjnIsQxm2nao6lOTwrQa7gKur6tYklzbLrwL2AhcA+4HvAD+10bbNrt8GvD/JG4CvAD/WV5aU3zMrSerhZShJUi8rC0lSr21RWRzJ7e6bWIYzkvxZktuS3JrkTR3rvCzJg0luaqZfnkM57kry+Wb/+zqWL+JcPGfsNd6U5KEkl7fWmcu5SHJ1kvuT3DI2b6qhCzbrGxrXKcO/SfLF5pxfm+TJ62y74e/vCMvwliR3j53zC9bZdtO+qXKdcvzBWBnuSnLTOttuyrnQglTVUk+MGma+BPwd4Bjgc8BZrXUuAD7EqF/xi4FPz6EcpwIvbB4/CfjLjnK8DPjgnM/HXcApGyyf+7no+P18DfieRZwL4IeAFwK3jM27EriieXwF8GtD3kdHWIZXA0c3j3+tqwzT/P6OsAxvAd48xe9rU87DeuVoLX8H8MvzPBdOi5m2Q7I4ktvdN01V3VtVn20eP8xoEMQ1X9q0BOZ+LlpeAXypqta7w35TVdXHgW+0Zk8zdME076PBZaiqj9RokEyATzHq0z4365yHaWzaeegrR5IAPw68d+j+tTy2Q2UxzTfrTfXte5slyTOBFwCf7lj8A0k+l+RDSZ43h8MX8JEkn0lyScfyhZ4LRn231/tnMO9zcdjE0AVA19AFizwvP80o3XXp+/0dqcuaS2FXr3M5bpHn4QeB+6rqjnWWz/tcaBNth8riSG5333RJTgD+CLi8qh5qLf4so8sxzwd+C/iPcyjCS6rqhYxGkvzZJD/ULmLHNvM6F8cAPwL8YcfiRZyLWSzkvCT5l8Ah4D3rrNL3+zsS7wKeDZwN3MvoEtCaInbMm1f/+deycaqY57nQJtsOlcWR3O6+qZLsZlRRvKeq/ri9vKoeqqpHmsd7gd1JNnUwsqq6p/l5P3Ato8sK4xb5TYTnA5+tqvs6yjn3czFmmqEL5n5eklwM/CPgJ6uq8x/wFL+/warqvqpaqapV4N3r7HtRfytHA/8E+IP11pnnudDm2w6VxZHc7r5pmuuvvw3cVlW/vs46T2/WI8m5jM7v1zexDMcnedLhx4waVW9prTb3czFm3U+O8z4XLdMMXTDN+2iwJOcBvwj8SFV9Z511pvn9HUkZxtumfnSdfc/1PIx5JfDFqjrQtXDe50JzsNUt7NNMjHr4/CWjXhz/spl3KXBp8ziMvuTjS8DngXPmUIaXMorrNwM3NdMFrXJcBtzKqIfJp4C/v8ll+DvNvj/XHGdLzkVznCcy+uf/X43Nm/u5YFQ53QscZPQp+Q3AdzH6Apc7mp8nN+t+N7B3o/fRJpZhP6O2gMPvjavaZVjv97eJZfi95nd+M6MK4NR5nof1ytHM/53D74WxdedyLpwWMznchySp13a4DCVJ2mJWFpKkXlYWkqReVhaSpF5WFpKkXlYWkqReVhaSpF5WFlpaSf6bZlC845o7fm9N8n3rrPsLzXcjfC7J2xZdVmmn86Y8LbUkvwIcBzwBOFBVv9qxzvnA/wq8sqq+k+TkqhoyfLekdVhZaKk14xfdCPwNoyFDVjrWeQejcYjevejySY8XXobSsjsZOIHRtxMet846YX7DbEvCykLLbw+jS0zvYfR1pV0+Avx0kifC6Du5F1Q26XHj6K0ugLSeJK8HDlXV7yfZBfxFkh+uqv80vl5VfTjJ2cC+JI8Be4FfWnyJpZ3LNgtJUi8vQ0mSenkZSttGkv+a0Rf8jHu0ql60FeWRHk+8DCVJ6uVlKElSLysLSVIvKwtJUi8rC0lSr/8fRuL4KsyrrHEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# an example of calculating kinetic energy\n", "ke = 0.5*(grid.interp((ds.uo)**2, 'X', **bd) + grid.interp((ds.vo)**2, 'Y', **bd))\n", "ke[0,0].plot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" } }, "nbformat": 4, "nbformat_minor": 4 }