Skip to content

Instantly share code, notes, and snippets.

@matt-long

matt-long/o2_heave_integral_analysis.ipynb

Last active July 2, 2019 08:10
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save matt-long/d1c34be6e111c3b8215057bd45225ed5 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save matt-long/d1c34be6e111c3b8215057bd45225ed5 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
o2_heave_integral_analysis.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Examine depth-integrated oxygen anomalies"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import os\n",
"\n",
"import xarray as xr\n",
"import numpy as np\n",
"\n",
"import esmlab"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spin up dask cluster"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/glade/work/mclong/miniconda3/envs/analysis/lib/python3.7/site-packages/dask_jobqueue/config.py:12: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
" defaults = yaml.load(f)\n",
"/glade/work/mclong/miniconda3/envs/analysis/lib/python3.7/site-packages/distributed/deploy/local.py:138: UserWarning: diagnostics_port has been deprecated. Please use `dashboard_address=` instead\n",
" \"diagnostics_port has been deprecated. \"\n",
"/glade/work/mclong/miniconda3/envs/analysis/lib/python3.7/site-packages/distributed/bokeh/core.py:74: UserWarning: \n",
"Port 8787 is already in use. \n",
"Perhaps you already have a cluster running?\n",
"Hosting the diagnostics dashboard on a random port instead.\n",
" warnings.warn(\"\\n\" + msg)\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "ea49fa653a0041aba2ff81cf273f7e64",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"VBox(children=(HTML(value='<h2>NCARCluster</h2>'), HBox(children=(HTML(value='\\n<div>\\n <style scoped>\\n .…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from ncar_jobqueue import NCARCluster\n",
"cluster = NCARCluster()\n",
"cluster.scale(4 * 9) # Ask for 4 x 9 workers\n",
"cluster"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table style=\"border: 2px solid white;\">\n",
"<tr>\n",
"<td style=\"vertical-align: top; border: 0px solid white\">\n",
"<h3>Client</h3>\n",
"<ul>\n",
" <li><b>Scheduler: </b>tcp://10.12.205.28:44026\n",
" <li><b>Dashboard: </b><a href='http://10.12.205.28/proxy/41161/status' target='_blank'>http://10.12.205.28/proxy/41161/status</a>\n",
"</ul>\n",
"</td>\n",
"<td style=\"vertical-align: top; border: 0px solid white\">\n",
"<h3>Cluster</h3>\n",
"<ul>\n",
" <li><b>Workers: </b>0</li>\n",
" <li><b>Cores: </b>0</li>\n",
" <li><b>Memory: </b>0 B</li>\n",
"</ul>\n",
"</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<Client: scheduler='tcp://10.12.205.28:44026' processes=0 cores=0>"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from dask.distributed import Client\n",
"client = Client(cluster) # Connect this local process to remote workers\n",
"client"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load monthly z-coordinate data"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"xr_open = dict(decode_times=False, decode_coords=False, chunks={'time': 48})\n",
"keep_vars = ['TLAT', 'TLONG', 'KMT', 'TAREA', 'z_t', 'z_w', 'z_w_bot', 'dz', \n",
" 'time_bound', 'time', 'THICKNESS', 'Z', 'sigma']"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (d2: 2, nlat: 384, nlon: 320, time: 816, z_t: 60, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * z_t (z_t) float32 500.0 1500.0 2500.0 ... 512502.8 537500.0\n",
" * z_w (z_w) float32 0.0 1000.0 2000.0 ... 475012.0 500004.7 525000.94\n",
" * z_w_bot (z_w_bot) float32 1000.0 2000.0 3000.0 ... 525000.94 549999.06\n",
" * time (time) float64 9.092e+04 9.094e+04 ... 1.157e+05 1.157e+05\n",
"Dimensions without coordinates: d2, nlat, nlon\n",
"Data variables:\n",
" dz (z_t) float32 dask.array<shape=(60,), chunksize=(60,)>\n",
" TLAT (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" KMT (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" TAREA (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" TLONG (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" time_bound (time, d2) float64 dask.array<shape=(816, 2), chunksize=(48, 2)>\n",
" O2 (time, z_t, nlat, nlon) float32 dask.array<shape=(816, 60, 384, 320), chunksize=(48, 60, 384, 320)>\n",
" PD (time, z_t, nlat, nlon) float32 dask.array<shape=(816, 60, 384, 320), chunksize=(48, 60, 384, 320)>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"droot = '/glade/p/cesm/community/CESM-DPLE/CESM-DPLE_POPCICEhindcast'\n",
"case = 'g.e11_LENS.GECOIAF.T62_g16.009'\n",
"datestr = '024901-031612'\n",
"stream = 'pop.h'\n",
"\n",
"variable_list = ['O2', 'PD']\n",
"\n",
"ds_list = []\n",
"for variable in variable_list:\n",
" dsv = xr.open_dataset(f'{droot}/{case}.{stream}.{variable}.{datestr}.nc', **xr_open)\n",
" ds_list.append(dsv)\n",
"\n",
"dsz = xr.merge(ds_list)\n",
"dsz = dsz.drop([v for v in dsz.variables if v not in variable_list+keep_vars])\n",
"dsz.attrs = {}\n",
"dsz"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load monthly sigma-coordinate data"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (d2: 2, nlat: 384, nlon: 320, sigma: 76, time: 816, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * sigma (sigma) float64 23.4 23.45 23.5 23.55 ... 27.0 27.05 27.1 27.15\n",
" * time (time) float64 9.092e+04 9.094e+04 ... 1.157e+05 1.157e+05\n",
" * z_w (z_w) float32 0.0 1000.0 2000.0 ... 475012.0 500004.7 525000.94\n",
" * z_w_bot (z_w_bot) float32 1000.0 2000.0 3000.0 ... 525000.94 549999.06\n",
"Dimensions without coordinates: d2, nlat, nlon\n",
"Data variables:\n",
" TLAT (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" KMT (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" TAREA (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" TLONG (nlat, nlon) float64 dask.array<shape=(384, 320), chunksize=(384, 320)>\n",
" time_bound (time, d2) float64 dask.array<shape=(816, 2), chunksize=(48, 2)>\n",
" Z (time, sigma, nlat, nlon) float32 dask.array<shape=(816, 76, 384, 320), chunksize=(48, 76, 384, 320)>\n",
" THICKNESS (time, sigma, nlat, nlon) float32 dask.array<shape=(816, 76, 384, 320), chunksize=(48, 76, 384, 320)>\n",
" O2 (time, sigma, nlat, nlon) float32 dask.array<shape=(816, 76, 384, 320), chunksize=(48, 76, 384, 320)>"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"droot = '/glade/p/cgd/oce/projects/DPLE_O2/sigma_coord/CESM-DPLE_POPCICEhindcast'\n",
"stream = 'pop.h.sigma'\n",
"\n",
"variable_list = ['O2']\n",
"ds_list = []\n",
"for variable in variable_list:\n",
" dsv = xr.open_dataset(f'{droot}/{case}.{stream}.{variable}.{datestr}.nc', **xr_open)\n",
" ds_list.append(dsv)\n",
"\n",
"dss = xr.merge(ds_list)\n",
"dss = dss.drop([v for v in dss.variables if v not in variable_list+keep_vars])\n",
"\n",
"dss['Z'] = dss.Z.where(dss.O2.notnull())\n",
"dss['THICKNESS'] = dss.THICKNESS.where(dss.O2.notnull())\n",
"\n",
"\n",
"dsigma = dss.sigma.diff('sigma').values\n",
"np.testing.assert_allclose(dsigma[0], dsigma[1:])\n",
"dsigma = dsigma[0]\n",
"sigma_bounds = np.empty((len(dss.sigma), 2))\n",
"sigma_bounds[:, 0] = dss.sigma.values-dsigma/2.\n",
"sigma_bounds[:, 1] = dss.sigma.values+dsigma/2.\n",
"\n",
"\n",
"dss.attrs = {}\n",
"dss "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## TESTING\n",
"\n",
"**subset domain for testing purposes**"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"dsz = dsz.isel(nlat=slice(200, 210), nlon=slice(200, 210)).compute()\n",
"dss = dss.isel(nlat=slice(200, 210), nlon=slice(200, 210)).compute()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute annual means"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (d2: 2, nlat: 10, nlon: 10, sigma: 76, time: 68, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * sigma (sigma) float64 23.4 23.45 23.5 23.55 ... 27.0 27.05 27.1 27.15\n",
" * time (time) float64 9.107e+04 9.143e+04 ... 1.152e+05 1.155e+05\n",
" * z_w (z_w) float32 0.0 1000.0 2000.0 ... 475012.0 500004.7 525000.94\n",
" * z_w_bot (z_w_bot) float32 1000.0 2000.0 3000.0 ... 525000.94 549999.06\n",
"Dimensions without coordinates: d2, nlat, nlon\n",
"Data variables:\n",
" time_bound (d2, time) float64 9.088e+04 9.125e+04 ... 1.153e+05 1.157e+05\n",
" Z (time, sigma, nlat, nlon) float64 1.165e+04 ... 7.388e+04\n",
" THICKNESS (time, sigma, nlat, nlon) float64 109.7 109.9 ... 5.025e+03\n",
" O2 (time, sigma, nlat, nlon) float64 147.3 147.4 ... 10.76 10.23\n",
" TLAT (nlat, nlon) float64 3.624 3.624 3.623 ... 6.052 6.051 6.05\n",
" KMT (nlat, nlon) float64 58.0 59.0 59.0 60.0 ... 59.0 57.0 56.0 56.0\n",
" TAREA (nlat, nlon) float64 3.743e+13 3.743e+13 ... 3.745e+13 3.743e+13\n",
" TLONG (nlat, nlon) float64 185.6 186.7 187.8 ... 193.4 194.6 195.7\n",
"Attributes:\n",
" history: \\n2019-06-21 16:28:16.231517 esmlab.resample(<DATASET>, freq=\"a..."
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dss = esmlab.resample(dss, freq='ann')\n",
"dss"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/glade/work/mclong/miniconda3/envs/analysis/lib/python3.7/site-packages/xarray/core/nanops.py:159: RuntimeWarning: Mean of empty slice\n",
" return np.nanmean(a, axis=axis, dtype=dtype)\n"
]
},
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (d2: 2, nlat: 10, nlon: 10, time: 68, z_t: 60, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * z_t (z_t) float64 500.0 1.5e+03 2.5e+03 ... 5.125e+05 5.375e+05\n",
" * time (time) float64 9.107e+04 9.143e+04 ... 1.152e+05 1.155e+05\n",
" * z_w (z_w) float32 0.0 1000.0 2000.0 ... 475012.0 500004.7 525000.94\n",
" * z_w_bot (z_w_bot) float32 1000.0 2000.0 3000.0 ... 525000.94 549999.06\n",
"Dimensions without coordinates: d2, nlat, nlon\n",
"Data variables:\n",
" time_bound (d2, time) float64 9.088e+04 9.125e+04 ... 1.153e+05 1.157e+05\n",
" O2 (time, z_t, nlat, nlon) float64 199.7 199.9 199.9 ... nan nan\n",
" PD (time, z_t, nlat, nlon) float64 1.022 1.022 1.022 ... nan nan\n",
" dz (z_t) float32 1000.0 1000.0 1000.0 ... 24996.244 24998.11\n",
" TLAT (nlat, nlon) float64 3.624 3.624 3.623 ... 6.052 6.051 6.05\n",
" KMT (nlat, nlon) float64 58.0 59.0 59.0 60.0 ... 59.0 57.0 56.0 56.0\n",
" TAREA (nlat, nlon) float64 3.743e+13 3.743e+13 ... 3.745e+13 3.743e+13\n",
" TLONG (nlat, nlon) float64 185.6 186.7 187.8 ... 193.4 194.6 195.7\n",
"Attributes:\n",
" history: \\n2019-06-21 16:28:17.339507 esmlab.resample(<DATASET>, freq=\"a..."
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dsz = esmlab.resample(dsz, freq='ann')\n",
"dsz"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute climatologies"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (d2: 2, nlat: 10, nlon: 10, z_t: 60, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * z_t (z_t) float64 500.0 1.5e+03 2.5e+03 ... 5.125e+05 5.375e+05\n",
" * z_w (z_w) float32 0.0 1000.0 2000.0 ... 475012.0 500004.7 525000.94\n",
" * z_w_bot (z_w_bot) float32 1000.0 2000.0 3000.0 ... 525000.94 549999.06\n",
"Dimensions without coordinates: d2, nlat, nlon\n",
"Data variables:\n",
" time_bound (d2) float64 1.031e+05 1.035e+05\n",
" O2 (z_t, nlat, nlon) float64 201.3 201.5 201.7 ... nan nan nan\n",
" PD (z_t, nlat, nlon) float64 1.022 1.022 1.022 ... nan nan nan\n",
" dz (z_t) float32 1000.0 1000.0 1000.0 ... 24996.244 24998.11\n",
" TLAT (nlat, nlon) float64 3.624 3.624 3.623 ... 6.052 6.051 6.05\n",
" KMT (nlat, nlon) float64 58.0 59.0 59.0 60.0 ... 59.0 57.0 56.0 56.0\n",
" TAREA (nlat, nlon) float64 3.743e+13 3.743e+13 ... 3.745e+13 3.743e+13\n",
" TLONG (nlat, nlon) float64 185.6 186.7 187.8 ... 193.4 194.6 195.7\n",
"Attributes:\n",
" history: \\n2019-06-21 16:28:17.339507 esmlab.resample(<DATASET>, freq=\"a..."
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dsz_clim = dsz.mean('time')\n",
"dsz_clim"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (d2: 2, nlat: 10, nlon: 10, sigma: 76, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * sigma (sigma) float64 23.4 23.45 23.5 23.55 ... 27.0 27.05 27.1 27.15\n",
" * z_w (z_w) float32 0.0 1000.0 2000.0 ... 475012.0 500004.7 525000.94\n",
" * z_w_bot (z_w_bot) float32 1000.0 2000.0 3000.0 ... 525000.94 549999.06\n",
"Dimensions without coordinates: d2, nlat, nlon\n",
"Data variables:\n",
" time_bound (d2) float64 1.031e+05 1.035e+05\n",
" Z (sigma, nlat, nlon) float64 1.217e+04 1.215e+04 ... 7.283e+04\n",
" THICKNESS (sigma, nlat, nlon) float64 93.88 93.62 ... 5.105e+03 5.081e+03\n",
" O2 (sigma, nlat, nlon) float64 149.2 149.5 149.8 ... 11.87 11.1\n",
" TLAT (nlat, nlon) float64 3.624 3.624 3.623 ... 6.052 6.051 6.05\n",
" KMT (nlat, nlon) float64 58.0 59.0 59.0 60.0 ... 59.0 57.0 56.0 56.0\n",
" TAREA (nlat, nlon) float64 3.743e+13 3.743e+13 ... 3.745e+13 3.743e+13\n",
" TLONG (nlat, nlon) float64 185.6 186.7 187.8 ... 193.4 194.6 195.7\n",
"Attributes:\n",
" history: \\n2019-06-21 16:28:16.231517 esmlab.resample(<DATASET>, freq=\"a..."
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dss_clim = dss.mean('time')\n",
"dss_clim"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Find `sigma_z_clim`: the isopycnal surface at `z = h` in the climatology\n",
"\n",
"We are using `z_w_bot`, with is the depth coordinate at the base of the tracer grid cell."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"k = 35\n",
"h = 671 m\n"
]
}
],
"source": [
"h_target = 630\n",
"\n",
"# find logical index at depth just below target\n",
"k = np.where(dsz.z_t*1e-2 >= h_target)[0][0]\n",
"\n",
"# get that depth\n",
"h = dsz.z_w_bot[k]\n",
"print(f'k = {k}')\n",
"print(f'h = {h.values/100:0.0f} m')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the `sigma_z_clim` as the mean of the layers above and below `z_w_bot[k]`"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2ae066662668>"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEGCAYAAABCa2PoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3debxcVZ3v/c83JwPzIGGQkO6gMqqMaUARuw1DQ6ShG3ggYCOot7lPX1FQTGPz6mu3eLmOCNgidhjjlQvYgoo+UciDiIKIJBEJmAi5jAEEoyjIkOGc7/1jrQo7lao6+5xU1Tmn9u/9eu1XTu1aa6+9K0n9zpplmxBCCGFDjRvpGwghhNAbIqCEEEJoiwgoIYQQ2iICSgghhLaIgBJCCKEtxo/0DRT1bbqpx7/udR0vRwMdL+K1sro1iK6/S+UA6l5R3dXNAY9dKqtr//66XNZLv1++wva2G3KNv37Xpv7d78v9x1l4/8pbbB+5IeVVwagKKONf9zqmfPQjnS/n5e59Jfa92p1yxr/UnXIAxnUxIHfzS76bv2iMW9Wdcsa/2r0PcNya7pV1z//+2OMbeo3f/b6fn9/yZ6XS9r3+4ckbWl4VjKqAEkII3WJggG7+dtT7IqCEECrJmNXuYltxBURACSFUVtRQ2isCSgihkozpj6Wn2ioCSgihsga6OrSv90VACSFUkoH+CChtFQElhFBZUUNpr5gpH0KoJAOr7VLHYCRNlXS7pCWSHpR0Vj5/g6T78vGYpPua5L9K0nOSHqg7/zpJ8yU9nP/cOp+XpC9JWibpfkn7bfgnsuEioIQQKsmY/pJHCWuAc2zvARwEfFDSnrZPsr2P7X2AG4GbmuS/Bmg0E//jwG22dwFuy68BjgJ2yccZwGUlH7ujIqCEEKrJ0F/yGPRS9jO2F+WfXwSWAFNq70sScCJwXZP8PwZ+3+CtY4G5+ee5wN8Wzn/Nyc+ArSS9vsxjd1IElBBCJaWZ8uUOYLKkBYXjjGbXlTQN2Be4p3D6EOBZ2w8P8Ta3t/0MpKAFbJfPTwGeLKRbTiGAjZTolA8hVJToL7/U6Qrb0we9orQZqWnrbNsvFN46mSa1k2FqdOMjPsIgAkoIoZJSp3z7FoqVNIEUTK61fVPh/HjgOGD/YVz2WUmvt/1MbtJ6Lp9fDkwtpNsJeHp4d94+0eQVQqikNA9FpY7B5D6SK4Eltr9Y9/ZhwFLby4dxmzcDp+WfTwO+Uzj/3jza6yDgj7WmsZEUASWEUFkDVqmjhIOBU4EZhWHCM/N7s6hr7pK0o6R5hdfXAXcDu0laLukD+a3PAIdLehg4PL8GmAc8AiwDLgf+27A+gDaLJq8QQiXVaihtuZZ9J032nrN9eoNzTwMzC69PbpL3d8ChDc4b+OAwb7djIqCEECrJiP5opGmrCCghhMoq2ZwVSoqAEkKoJCNWuW+kb6OnREAJIVRSmtgYTV7tFAElhFBZ7eqUD8moCigTXobtFnR+suf4V7u37ef4l7qzZ3Xfq93bG1sD3ZuQ63Hd+w+vLs4z1uru/Bsct2pNV8oByi16NYrYot9RQ2mnURVQQgihmwaihtJWEVBCCJWUOuXjK7Cd4tMMIVRSdMq3XwSUEEJl9cc8lLaKgBJCqKSYKd9+EVBCCJU1EKO82qqjn6akj0h6UNIDkq6TtFEnywshhLLS4pDjSh2hnI59UpKmAB8Gptt+C9BHWsY5hBBGnBGr3VfqCOV0uslrPLCxpNXAJoyCHcVCCAHAJiY2tlnHPk3bTwFfAJ4AniHtKHZrfTpJZ0haIGnB6ldf6tTthBBCHTFQ8gjldLLJa2vgWGBnYEdgU0l/X5/O9hzb021Pn7DRpp26nRBCWIdJNZQyx2AkTZV0u6Qlud/4rHz+hsIOjo9Juq9J/iMl/VrSMkkfL5yfIWlR7oeem/enR9LWkr4l6X5JP5f0lvZ8Khumk/W9w4BHbf/W9mrgJuDtHSwvhBCGpI2d8muAc2zvARwEfFDSnrZPsr2P7X2AG0nfg+uQ1AdcChwF7AmcLGlPSeOAucCs3A/9OK/tL38ecJ/tvYD3Apds0AfRJp0MKE8AB0naRJJI21gu6WB5IYRQmim3n3yZTbhsP2N7Uf75RdJ33ZTa+/k78ETq9pbPDgCW2X7E9irgelLrzjbAStsP5XTzgePzz3sCt+XylgLTJG0/9E+hvTrWKW/7HknfBBaRovcvgDmdKi+EEIbCwOoOrOUlaRqwL3BP4fQhwLO2H26QZQrwZOH1cuBAYAUwQdJ02wuAE4CpOc0vgeOAOyUdAPw5sBPw7Abc934lkq22vbjZmx0d5WX7X4F/7WQZIYQwPBrKfiiTJS0ovJ5je71fkCVtRmraOtv2C4W3TqZx7STdyPps25JmARdJmgTcSvrlHOAzwCW5T2Yx6Rf2Dd2r4A7g3ib3U7MzMK3ZmzFTPoRQSWZIM+VX2J7eKoGkCaRgcq3tmwrnx5NqE/s3ybqc12oekGoaTwPYvptUu0HSEcCu+fwLwPvyeQGP5mND3Gt7RqsEkn7Y6v0YhB1CqKz+XEsZ7BhM/lK/Elhi+4t1bx8GLLW9vEn2e4FdJO0saSJpAvjN+brb5T8nAecCX82vt8ppAf4L8OO6GtGQDRZMyqSJgBJCqCRbDHhcqaOEg4FTgRmFYcIz83uzqGvukrSjpHnpPrwGOBO4hdSZ/w3bD+aksyUtAe4Hvmu7VkPYA3hQ0lLS6LCzhv1BrHtfO0jaIf+8raTjJL25bP5o8gohVFLqlG/Psiq276RJ34Pt0xucexqYWXg9D5jXIN1sYHaD83cDuwz/jtcn6b8CH08/6rPA6cCDwKclfc72lYNdIwJKCKGiYk/5OmcCbwY2Js15eZPt3+RJ6reTmvRaGlUBpe+VNWy1+PmOl6NXVnW8jLVefqU75azZ0AEeQzDg7pU1vosL843r4pfLwEB3ynEX/67GmNQpH8uqFKy2/TLwsqT/Y/s3ALafl1TqH9KoCighhNBNsTT9OgYkTcgrm7y7djJvO1Lqg4qAEkKopNpM+bDWcaSKG3Uj0rYBzilzgQgoIYTKGogaylq2n2hy/ingqTLXiE8zhFBJNqweGFfqqBJJpw83b7U+qRBCyFKTV9vmofQESZ8gTcQclmjyCiFU1hDW8up5kuYAmwPvGe41IqCEECophg2v5xTgANvDHtMeASWEUFGqTHOWpGOAd+aXd9j+boNkRwPfkHSs7f8znHKq8WmGEEIDVdhTXtKnSWt9/SofH87n1mH7R6R1x74+3LKihhJCqKQ0yquLKzGMnHcD+9SasiTNJe2f8s/1CW0/IOn4+vNlRQ0lhFBJ7dwCeAzYqvDzls0S5f3tTxpuIVFDCSFU1lhvzirp08AvJN1OWhH5ncB5jRLa7pd0LHDRcAqKgBJCqKSqjPKyfZ2kHwF/QQoo59YWfmziLklfBm4AXipcZ9FgZUVACSFUVhVGeUm6zfah5F0g68418vb85/mFcwYG3dExAkoIoZJssaaHA0peJXgTYHLe06RWHdsC2LFZPtvvGm6ZvftphhDCINrVKS9pqqTbJS2R9KCks/L5GwpbAj8m6b4m+Y+U9GtJyyR9vHB+hqRFkh6QNFfS+Hx+S0nflfTLXN77Glz2vwILgd3zn7XjO8ClLZ5le0lXSvp+fr2npA8M+iEQASWEUFG1PpQ2jfJaA5xjew/gIOCDkva0fZLtfWzvA9wI3FSfMY+supS0N/yewMn5S3wcMBeYZfstpF0UT8vZPgj8yvbewF8BF0qauM7z2ZfY3hn4mO032N45H3vb/nKLZ7mGtL99rRbzEHB2mQ8hAkoIobLaFVBsP1PrtLb9IrAEmFJ7X5KAE4HrGmQ/AFhm+xHbq4DrgWNJ+5CstP1QTjcfqM0RMbB5vu5mwO9JQa3Rvf37oA+wrsm2vwEM5PxrgP4yGSOghBAqaYjzUCZLWlA4zmh2XUnTgH2BewqnDwGetf1wgyxTgCcLr5fncyuACZKm5/MnAFPzz18G9gCeBhYDZ23IGlx1XpK0DXmzLUkHAX8skzE65UMIlTWEeSgrbE8fLJGkzUhNW2fbfqHw1sk0rp0ADW/Cti1pFnCRpEnArbxWC/lr4D7SyKs3AvMl/aRYpqSDbd8laZLtlYPde8FHSSPC3ijpLmBb4P8pkzECSgihkmxY08bNsyRNIAWTa23fVDg/nrS97v5Nsi7ntZoHwE6kmge27ybVbpB0BLBrTvM+4DO2DSyT9Cip8/3nhet8KZd5N7DfEB7lQeAvgd1Iwe7XjMk95Veuhsef7ngxA6tWd7yM18pa1Z2C2lbbLUHdaykdN6GL/0S7+Fz0dacsjeviM40fXV8nZbRrYmPuy7gSWGL7i3VvHwYsrdunveheYBdJO5O22p1FWkoeSdvZfi7XUM4FLsh5ngAOBX4iaXvSl/8jddddLelqYIqkL9UXavvDTe7nbtv7kQJL7fkWUSIojb1/ASGE0Aa1PpQ2ORg4FVhcGBp8nu15pACxTnOXpB2BK2zPtL1G0pmkkVV9wFW2a1/msyUdTaohXGb7h/n8p4BrJC3mtdnvK+ru6WhSMJtBGi7ckqQdSH03G0val3XnrWxS5kOIgBJCqCy3KaDYvpPGfSHYPr3BuaeBmYXX84B5DdLNBmY3yX/EIPe0Arhe0hLbvxzkESD1y5xOanIr1rJeoMnaX/UioIQQKqsii0O+Iuk2YHvbb5G0F3CM7f9RTGR7LjBX0vG2bxxOQTFsOIRQSXZbJzaOZpeT9j5ZDWD7flIzXDN3xUz5EEIYEtE/MK7UMcZtYvvndecaToLMriZmyocQwtDYKnWMcSskvZHXJiqeADzTIv2wZ8p3tA9F0lbAFcBbSA/z/jyuOoQQRlRV9kMhrfs1B9hd0lPAo8Dft0g/amfKXwL8wPYJeeGyUkPPQgih45z6UXqd7UeAwyRtCozLa4210mim/AllyupYQJG0BWmrydMB8qJnXZrlF0IIg+v1UV55JeOtba+w/ZKkiZL+AfhoXhl5PbYXSVpnprztUrPBO1lDeQPwW+BqSXuTJtacZfulYqK8yNoZABtp0w7eTgghvMa5U75X5XXA/oPUhPUw8G/A/yLNzH9Pi3x9pDky00gx4ghJNFgBYD2dDCjjSVP1P2T7HkmXAB8H/nsxke05pPY9tuybXIEKaAhhtOjxJq9/Afa3vUzSfqQ1vWbZ/tYg+b4LvEpaxXhIazp1MqAsB5bbri3h/E1SQAkhhFGhB0ZwtbLK9jJY24z1aIlgArCT7b2GU2DHAort30h6UtJutn9NWsjsV50qL4QQhsLu+YCynaSPFl5vVnzdognr+5KOsH3rUAvs9CivDwHX5hFej5CWXA4hhFGhx4cNXw5s3uJ1Mz8DvpW3IF5N6pi37S0Gy9jRgGL7PmDQTWlCCGEk9HIfiu1PDjPrhcDbgMV5v5XSYnHIEEIlGTHQw6O8NsDDwANDDSYQASWEUGE9XEHZEM8AP8qLQ67dOnikhw2HEMLo1fud8kCaV2K71Fpc2aP5mJiP0iKghBCqq01VFElTga8BO5DmbsyxfYmkG0gzzgG2Av5ge58G+Y8kLVXVR9rJ8TP5/AzgC6Qv9oXAB/IOj7N5bXLieGAPYFvbv29we8skfRO42vagI203oO8lAkoIobraWENZA5yT53tsDiyUNN/2SbUEki6kwSKLeWb6pcDhpPl790q6GVgKzAUOtf2QpPOB04ArbX8e+HzO/zfAR5oEE4C9SPufXJFHbl0FXG/7hbr7uNj22ZK+S4NQa/uYwT6ECCghhEoyMDDQti2AnyEvCW/7RUlLSPuz/wpAkoATSfu71zsAWJYXcUTS9cCxpKWrVtp+KKebT9oo68q6/CdTt2d93b29SBoyfLmkd+a0F+Vay6dqkx9Jy7JAqhENy+gKKDZeuXLwdBtaTP9QmhM3jPr6ulNQF9uCNXFIzaobVtb4Lv4TndC9srr2XF38u2JSF8t6vg3XMEP5fzNZ0oLC6zl52aj1SJoG7AvcUzh9CPCs7YcbZJkCPFl4vRw4EFgBTJA03fYC0oq/U+vK2gQ4Ejiz2Y3nGtC7SfMAp5GGBV+b72kesCuA7YU5yz62L6m7xlnAHc3KqBldASWEELpoCANjV9gedE6dpM2AG4Gz65qUWtUiGkU123Ze4PEiSZOAW1l/p8W/Ae5q0dwFaRjw7cDnbf+0cP6bucZS7zRSf07R6Q3OrScCSgihuto4bljSBFIwudb2TYXz44HjgP2bZF3OujWPnYCnAfKGhIfk6xxBrk0UzKJFc1f2Xtt31t3rwbbvsv3hwrmTgVOAnXMfTs3mwO8GKQOIgBJCqKz2be+b+0iuBJY0mK9xGLDU9vIm2e8FdpG0M/AUKUickq+7ne3ncg3lXOCCQplbAn9J690XAb5EWvm96N8bnPspqR9oMqlZrOZF4P5BygAioIQQqqx9NZSDgVOBxZLuy+fOsz2PBrUISTuShgfPzMOAzwRuIQ0bvsr2gznpbElHA+OAy2z/sHCZvwNurd9jqlDG24C3A9vWLRK5RS5nHbYfBx4nLbsyLBFQQgjVZHD7RnndSeO+EGyf3uDc06RNrGqv55E6yOvTzQZmN7nuNcA1LW5rIrAZ6Xu+uCjkC7TY0lfSccBnge1IzzQ6FocMIYTRrXdnytu+A7hD0jW59lHW54C/sb1kqGVGQAkhVFcPL+ZVm6gIfFnSUCYqPjucYAIRUEIIVdbDAYXhT1RckJeM+TbrLg55U/MsSQSUEEI1DW1i45hTm6iYm74AkLQ1MNV2q1FbWwAvA0cULwdEQAkhhGZ6eYOtGkk/Ao4hfd/fB/xW0h22P9oove1h76wbu8uEEKprQOWOsW3LPGv/ONKKw/uT5sY0JGlXSbdJeiC/3kvSv5QpKAJKCKGy5HLHGDde0utJi1N+r0T6y0mLUK4GyM1js8oUFAElhFBNHsIxtp1PmjS5zPa9kt5AWt+rmU1s/7zuXP0aYg2V6kORNMn2ysHOhRDC2KGe7pSvsf2fwH8WXj8CHN8iywpJbySHUkknkJfmH0zZTvm7WX/dl0bnQghh7Bj7tY9BSdoW+AfS0vVrv/Ntv79Jlg8Cc4DdJT1F2g54sPXCgEECiqQdSGv1byxpX16bVroFsEmZAkIIYdQaGOkb6IrvAD8B/n9g0M2gcg3mMEmbAuPyBl2lDFZD+WvSOvg7AcUVNF8EzitbSAghjDo9Pg+lYBPb55ZNLOl/Ap+z/Yf8emvS9saDjvRq2Slve67tdwGn235X4TimzKzJEEIYzSoyyut7kmYOnmyto2rBBMD28xQWsmylVB+K7RslvRt4M7BR4fz5Q7jJEEIYXcZ+sCjjLOA8SauAVQy+enBfcdCVpI2BSWUKKjvK66ukPpN3AVeQlj6uH1YWQghhlLG9+eCp1vF14DZJV5NC7vuBuWUylh3l9Xbbe0m63/YnJV1IiXVdQghhNOuB5qxB5d0k3wPsbPtTkqYCr28w1wQA25+TtBg4lFSb+ZTtW8qUVTagvJL/fDnvNPY7YOeSecsbP55xk7dp+2XXM66L8zn7ulRW33obsPVEWZ7UveXmPLF7ZfV3qayBSd37u+qf1MX/V8vacA3TtmVV8pf014AdSGPH5ti+JK/au1tOthXwB9v7NMh/JHAJaSfFK2x/Jp+fQVoteCKwEPiA7TX5vb8CLgYmACts/2WT2/tKvqcZwKeAPwGXAn9Rdw+3AD8Avm/7+8D3h/o5lP1X/T1JWwGfBxaR/iquGGphIYQwqrSvhrKGNBJqkaTNgYWS5ts+qZYgt+z8sT6jpD7SF/zhwHLgXkk3A0tJTU2H2n5I0vnAacCV+fv4K8CRtp+QtF2LezvQ9n6SfgGpk13SxAbpTgOOBP5N0q7APaQAc5vtP5X5EMp2yn8q/3ijpO8BG9le74MJIYSxpF1NXrafIc8mt/2ipCWkOXy/grXNTieSagn1DiAti/JITns9cCzwW2Cl7YdyuvmkNbauBE4BbrL9RC7zuRa3tzoHrdrM921pMAPH9m9IWwpfI2kccCBwFPBPkl4h7V//uVafw2ATG49r8V6pDVdCCGHUKh9QJktaUHg9x/acRgklTQP2Jf2GX3MIaSfERmtoTQGeLLxeTvoyXwFMkDTd9gLSYKipOc2u+b0fkfaLv8T215rc+5eAbwHbSbogX6flnBLbA6TVUO4GPiFpMmleYkuD1VD+pr6c/KcoueFKCCGMWuUDygrb0wdLJGkz4Ebg7LxkfM3JwHXNsjW6M9uWNAu4SNIk4FZeW6RxPLA/qeN8Y+BuST8r1GaKF7pW0kJe62T/21Zb/EraGfgw8Oesu1RLsy2D12oZUGobrUjaiLSY2LRCngqMjwgh9Kp2T1qUNIEUTK4ttt5IGk/ai2T/JlmX81rNA9LKJE8D2L6bVLtB0hGkmkktzwrbLwEvSfoxsDewXkCR9FZgd+A5YEmJ/eK/TWpWu5khLk5TtlP+28AfSB3yr+ZzEVBCCGNb+0Z5ifQlvMT2F+vePgxYant5k+z3ArvkmsFTpL1HTsnX3c72c7mGci5wQc7zHeDLOVhNJDWRXVR3T1vmdFOB+0m1k7dKegI4tq4GVfSq7S+VfPR1lA0oO9k+cjgF5M6gBcBTto8ezjVCCKET2lhDORg4FVgs6b587jzb80gBYp3mrjz94grbM22vkXQmac+SPuAq2w/mpLMlHU1aJusy2z8EsL1E0g9IgWIgX+uBunv6FOm7d0buE6l9H3+aFJg+1ORZLpH0r6QmtrVblNheNNiHUDag/FTSW20vLpm+6CxgCWmF4hBCGD3aN8rrThr3hWD79AbnnqawPlYOPPMapJsNzG5y3c+TpnI0cxiwVy2Y5Dz9ks4DWn2Xv5UUHGfwWpOXaTxCbR1lA8o7gNMlPUqKWLW1YPZqlUnSTsC7SdHwoyXLCiGEzuuNhR9bWVWbBFmUa0StNkf8O+ANtlcNtcCyAeWooV44uxj4J9KwtoYknQGcAbBR31CXnAkhhA3Q2wFlo7p9rGpE68Uef0ma1d9qbktDZSc2Pj7UC+d2v+dsL8xLBDS79hzS7mBsOXH73v7rDSGMKurtDbaeYd19rIp+0yLf9sBSSfeybh/Khg0b3kAHA8fkdfg3AraQ9HXbpbaSDCGEMHx5L6vh+NfhltmxgGL7n0nLBNQWMftYBJMQwqgSbSLrsX3HcPN2cXnQEEIYRUru1tjjHffrkXSQpHsl/UnSKkn9kprNWVlHV9bQtv0j4EfdKCuEEEqrWLAo6cukuTP/CUwH3gvsUiZj9zaACCGE0aaHA4qk/Vq932qiou1lkvps9wNXS/ppmTIjoIQQKkn0/CivC1u812qi4st5v5T7JH2ONFps0zIFRkAJIVRTj/ePbMAor1NJ/etnAh8hrQV2fJmMEVBCCNXVwwGlSNJbgD1JUzgAaLR/Sl7r64I8IvdV4JNDKScCSgihuioQUPJCj39FCijzSCuf3AmsF1DyWl/bSprYyaVXusIbTWDlbjt2vpzx7VmyuoyBvu6U1dVnmtC9svonda+sNV0sq39id8rqb7XARpsNNNqlvFNuac9lernJq+AE0l4pv7D9PknbA1e0SP8YcFfe1/6l2skGy/KvZ1QFlBBC6KpqBJRXbA9IWiNpC9IaXW9okf7pfIyjxTqMjURACSFUk3t+lFfNAklbAZcDC4E/AT9vltj2kPpNiiKghBCqqwI1FNv/Lf/41bwp1xa272+WXtJ3Wf+T+SNps67/sP3q+rmSWHolhFBZvbz0iqTd85/71Q7gdcD4QSY9PkKqxVyejxeAZ0n72V/eqsyooYQQqqtNwULSVNKoqR1IuxzOsX2JpBuA3XKyrYA/2N6nQf4jgUtIWwBfYfsz+fwM4AukfeMXAh/IG2T9FWm/+EfzJW6yfX7dZT9K2muq0QTHVhMb97X9zsLr70r6se13SnqwSR4gAkoIoapMO5u81gDn2F4kaXNgoaT5tk+qJZB0IanpaB157selwOHAcuDePMJqKTAXONT2Q5LOB04DrsxZf2L76GY3ZPuM/ONR9c1UkjZqkKVmW0l/ZvuJnPbPgMn5vZZDiaPJK4RQSaJ9TV62n6mtjWX7RWAJMGVtWZKAE4HrGmQ/AFhm+5E89+N64FhgG2Cl7YdyuvmUnLFep9E6XK3W5joHuFPS7ZJ+BPwEmC1pU1KAaypqKCGEyhpC/8hkSQsKr+fk3WbXv6Y0DdgXuKdw+hDgWdsPN8gyBXiy8Ho5cCCwApggabrtBaT5JFML6d4m6ZekIb4fs71Oc5SkHfK1N67bCngLYJNmD2p7nqRdgN1znqWFGs7FzfJBBJQQQpWVDygrbE8fLJGkzYAbgbNtF/cQOZnGtRNYf893ANu2pFnARZImAbeSmtYAFgF/bvtPeVfcb7P+EvN/DZwO7MS6WwG/CJzX4N73K9SyVpL2lm+appEIKCGE6mrjCC5JE0jB5FrbNxXOjweOA/ZvknU569Y8diLVOrB9N6l2g6QjSCOtKAarXKP4iqTJtlcUzs8F5ko63vaNJR7h6tzZ32oZhytJta+GIqCEEKqpjUOCcx/JlcCSBkuUHEZqNlreJPu9wC6SdgaeIm1udUq+7na2n8s1lHOBC/L5HUhNaJZ0AKk//HdNrv89SacA0yh85zcYFbYlaSRZq4Dy2xbvRUAJIVRY+2ooB5OWfV8s6b587jzb80gBYp3mLkk7koYHz8zDgM8krVDWB1xV6A+ZLeloUsC4zPYP8/kTgH+UtAZ4BZhlu9nTfIc0umwhsLLZA9ieNqQnbiACSgihstq19IrtO2nym73t0xucexqYWXg9j7QScH262cDsBue/TNqqt4ydbB9ZMu0GiWHDIYTK6uWZ8gU/lfTWbhQUNZQQQjW1d2LjaPYO4HRJj5KavEQaRbZXuwuKgBJCqK5qBJSjhpI4DzB4D/AG2+fnmfI72G66QnFNNHmFECqpnTPlRzPbj5OGJc/IP79M6+/+rwBvI82dgTRv5dIyZUUNJYRQWRoY49GihLwF8HTSIpVXAxOAr5NGpjVyoO39JP0CwPbzkkrtxxk1lBBCNXkIx9j2d8Ax5O188wizVjsxrs4LVhpA0mYVRvAAABAOSURBVLakFZQHFQElhFBZVWjyAlblOSq1ALHpIOm/BHwL2E7SBcCdwP8sU1A0eYUQqmvsB4syviHpP4CtJP0D8H7gimaJbV8raSFwKKmr6W9tLylT0KgKKGs2Fiv2ntTxctzFetlAlz7hgQndKafbZfV3/p/DWgOTuvftMjCpvyvleFL3Nk3XhLG3QXsP1D4GZfsLkg4n7by4G/AJ2/ObpZd0EPCg7Uvz680lHWj7nmZ5aqLJK4RQXRXoQ5H0Wdvzbc+2/THb8yV9tkWWy0hbANe8lM8NKgJKCKGanJZeKXOMcYc3ONdqboqK64LZHqBka1YElBBCJfX6PBRJ/yhpMbCbpPvzsTjPmL+/RdZHJH1Y0oR8nAU8UqbMUdWHEkIIXdV0gd6e8L+B7wOfBj5eOP+i7d+3yPf/kkZ6/Qupwe824IwW6deKgBJCqKyxWvsow/YfScvWnyxpb/JGXaQ94psGFNvPkZbcH7IIKCGEauqBDvcyJH2YVMOo7SL5dUlzbP97k/QbAR8A3gxsVDtv+/2DldWxPhRJUyXdLmmJpAdzO1wIIYwaFemU/y+k5VQ+YfsTwEHAP7RI/7+AHUh70t9B2pL4xTIFdbJTfg1wju09SA/wQUl7drC8EEIYknYFlGa/QEu6QdJ9+XissJtjff4jJf1a0jJJHy+cnyFpkaQHJM3N+9MX8/2FpH5JJ7S6PaA48amf1tv8vsn2fwdeyvvSvxsotZ9Kx5q8bD8DPJN/flHSEmAK8KtOlRlCCKWZdnbK136BXiRpc2ChpPm2T6olkHQhqU9jHXndrEtJw3uXA/dKuhlYCswFDrX9kKTzgdNIe9fX8n2WtHVwK1cD90j6Vn79t7VrNLE6//kHSW8BfkPaj35QXRk2LGkasC+w3kxLSWdIWiBpwZpXXurG7YQQAtC+YcO2n7G9KP/8IlD7BTqVk/YYOZG6veWzA4Blth+xvQq4HjgW2AZYafuhnG4+cHwh34eAG4HnBrm3LwLvI3XEPw+8z/bFLbLMkbQ1aZTXzaRKQKuJkGt1vFNe0makhz7b9gv179ueA8wB2GT7qRXoIgshjBrlv3EmS1pQeD0nf3etp8kv0IcAz9p+uEGWKcCThdfLgQOBFcAESdNtLwBOIO1rgqQppFWEZwB/0eQ+NiINAX4TsBj4iu01zR5Q0lm2LwGW2H4e+DHwhmbpG+loQJE0gRRMrrV902DpQwihW2oTG0taYXv6oNds/gv0yTSundRupZ5tW9Is4CJJk4BbSU1rABcD59ruT5WfhuaSmq9+QpoZvwdwdovbfx9wCfDvwH4t0jXVsYCSq3hXkqLdFztVTgghDIvd1g22mv0CnTvSjwP2b5J1Obnmke0EPJ1u0XeT549IOgLYNaeZDlyfg8lkYKakNba/XbjOnrbfmvNeCQy2he8SSY+Rlq0vzqQvvQd9J2soBwOnAosLIxvOsz2vg2WGEEJ5bYong/wCfRiw1PbyJtnvBXaRtDPwFGlS4Sn5utvZfi7XUM4FLgCwvXOh7GuA79UFE3itcx3ba1rUZGppTpa0A6mT/5iWiZvo5CivO2k9NC2EEEZUG2fKt/oFehZ1zV2SdgSusD0zf9mfSfoi7wOusv1gTjpb0tGkAVSX2f7hEO5pb0m1ZjcBG+fXtRrHFg3y/BZYnPeeH7KYKR9CqCYDbWryavULtO3TG5x7GphZeD0PWK/1xvZsYPYgZa93/Xy+r1W+Jnn6JU2WNDGPOBuSCCghhOqKcaWNPA7clefCrJ3LUaYvPAJKCKGyenlxyA3wdD7GAZsPJWMElBBCZbVzlFevsP3J4eaNgBJCqKaKrDY8VJJup8EnY3vGYHkjoIQQKilNbIyI0sDHCj9vRFrupekM+6JRFVD6NzEv7D3kgQVDpvHdW4+6b0L/4InaYEKXygHYeGKpf1ttsdmklV0ra/OJ3Strq4mvdKWcLSd0pxyAzcZ3/v9uzYXtutDYX5q+7WwvrDt1l6Q7yuQdVQElhBC6KWoo65P0usLLcaQZ/juUyRsBJYRQTdGH0sxC0icjUlPXo6QdHAcVASWEUFHtXcurVxSXdRmqruyHEkIIo5Jd7qiAvPvjDoXX75X0HUlfqmsGayoCSgihmlyZPeXL+g9gFYCkdwKfAb5G2mWy4d4v9aLJK4RQXRWpfZTUZ/v3+eeTSJuI3QjcWFjwsqWooYQQqsslj2roy3u3ABwKFFc2LlX5iBpKCKGyNFCd9qwSrgPukLQCeIW00yOS3kRq9hpUBJQQQjWZmNhYYPsCSbcBrwdutde2B44DPlTmGhFQQgiVJBwTG+vY/lmDcw+VzR99KCGE6mrTsGFJUyXdLmmJpAclnZXP3yDpvnw81qxzW9KRkn4taZmkjxfOz5C0SNIDkubW+jgkHSvp/nzdBZLe0aZPZINEDSWEUF3tq6GsAc6xvUjS5sBCSfNtn1RLIOlCGvRFSOoDLgUOB5YD9+bNrZYCc4FDbT8k6XzgNNLe9bcBN9u2pL2AbwC7t+thhitqKCGEaqr1oZQ5BruU/YztRfnnF4ElwJTa+5IEnEjd3vLZAcAy24/kbXevB44FtgFWFpqc5pNW/sX2nwp9HJsySsaiRUAJIVSWBgZKHcDk3LRUO85oek1pGrAvcE/h9CHAs7YfbpBlCvBk4fXyfG4FMEHS9Hz+BGBqoZy/k7QU+P+A9w/12TshmrxCCBU1pGVVVtiePlgiSZsBNwJn236h8NbJNK6dQFqEcb2by81Zs4CLJE0CbqWwL4ntbwHfyrPaPwUcVu5ROicCSgihmkxbZ8pLmkAKJtfavqlwfjxwHGkZ+EaWU6h5ADuR9nTH9t2k2g2SjgB2rc9s+8eS3ihpsu0V7XiW4YomrxBCdbWpDyX3kVwJLLH9xbq3DwOW2l7eJPu9wC6SdpY0EZgF3Jyvu13+cxJwLvDV/PpNuUwk7QdMBH5X6pk7KAJKCKGyZJc6SjgYOBWYURgmPDO/N4u65i5JO0qaB2B7DXAmcAupM/8bth/MSWdLWgLcD3zXdm05lOOBB/Iw5EuBkwqd9CMmmrxCCNXVpu9g23fSuC8E26c3OPc0MLPweh4wr0G62cDsBuc/C3x2+HfcGRFQQgjVZEN/rL3SThFQQgjVNfKtRD1lVAWUrTd5mRP3WdDxcjbpW9XxMmq27HulK+Vs3qVyALbse7lrZW3T96fulTWue8/1ur41gydqRznjJnalHICNNalrZV3YrgtFQGmrURVQQgihawzEnvJtFQElhFBRBkcfSjtFQAkhVJOJTvk2i4ASQqiu6ENpqwgoIYTqioDSVhFQQggVNaTFIUMJHV16pdkuZCGEMOIMDAyUO0IpHQsohV3IjgL2BE6WtGenygshhCFr0xbAIelkk9faXcgAJNV2IftVB8sMIYSSYumVdutkk1ezXcjWIemM2i5orzy/soO3E0IIBQZ7oNQRyulkQGm4C9l6J+w5tqfbnr7x1t1buiGEEBhwuSOU0skmr6a7kIUQwqgQ/SNt1cmAsnYXMuAp0iYzp3SwvBBCKM+OEVxt1rEmr0F2IQshhJHXplFekqZKul3SEkkPSjorn7+hsIPjY3mHxUb5G06xkDRD0iJJD0iam/enR9J7JN2fj59K2rtNn8gG6ejExma7kIUQwsgz7u9v18XWAOfYXiRpc2ChpPm2T6olkHQh8Mf6jIUpFoeTugrulXQzsBSYCxxq+yFJ5wOnkfaufxT4S9vPSzoKmAMc2K6HGa7YUz6EUE215evb0Clv+xnbi/LPL5JaZdaOapUk4ETq9pbP1k6xsL0KqE2x2AZYafuhnG4+aS95bP/U9vP5/M9IfdQjLgJKCKG6PFDuGAJJ04B9gXsKpw8BnrX9cIMszaZYrAAmSJqez5/AugOdaj4AfH9IN9khsZZXCKGSDLj8kODJkorbyc6xPac+kaTNgBuBs22/UHjrZBrXTqDJFAvbljQLuEjSJOBWUtNasbx3kQLKO8o+SCdFQAkhVJOHtMHWCtvTWyWQNIEUTK61fVPh/HjgOGD/JlmbTrGwfTepdoOkI4BdC9fdC7gCOMr278o+SCdFk1cIobLc31/qGEzuI7kSWGL7i3VvHwYstb28Sfa1UywkTSRNsbg5X3e7/Ock4Fzgq/n1nwE3AacW+lhGnDyKJvZI+i3w+BCzTSa1NfaaXnyuXnwm6M3nGu3P9Oe2t92QC0j6Aek5y1hh+8gW13oH8BNgMVCr9pxne56ka4Cf2f5qIf2OwBW2Z+bXM4GLgT7gKtsX5POfB44m/fJ/me2L8/krSB30te/LNYPVoLphVAWU4ZC0YDR8kO3Wi8/Vi88EvflcvfhMofOiySuEEEJbREAJIYTQFr0QUNYbutcjevG5evGZoDefqxefKXTYmO9DCSGEMDr0Qg0lhBDCKBABJYQQQluM2YDSbLnnsazZEti9QFKfpF9I+t5I30u7SNpK0jclLc1/Z28b6XtqB0kfyf/+HpB0naSNRvqewtgwJgNKYbnno4A9gZMl7Tmyd9UWtSWw9wAOAj7YI88FcBZpBdZecgnwA9u7A3vTA88naQrwYWC67beQJtrNGtm7CmPFmAwoNF/ueUwbbAnssUrSTsC7SesO9QRJWwDvJC23ge1Vtv8wsnfVNuOBjfMaVJsQW3eHksZqQGm23HPPaLIE9lh1MfBPvLYkRS94A/Bb4OrclHeFpE1H+qY2lO2ngC8ATwDPAH+0fevI3lUYK8ZqQGm43HPX76JDWiyBPeZIOhp4zvbCkb6XNhsP7EdaX2lf4CVgzPflSdqaVNvfGdgR2FTS34/sXYWxYqwGlKbLPY91zZbAHsMOBo6R9BipaXKGpK+P7C21xXJgue1aDfKbpAAz1h0GPGr7t7ZXk1a0ffsI31MYI8ZqQGm63PNYNsgS2GOS7X+2vZPtaaS/px/aHvO/8dr+DfCkpN3yqUOBX43gLbXLE8BBkjbJ/x4PpQcGG4TuGJMbbNleI+lM4BZeW+75wRG+rXY4GDgVWCzpvnzuPNvzRvCeQnMfAq7Nv9Q8ArxvhO9ng9m+R9I3gUWkUYe/IJZhCSXF0ishhBDaYqw2eYUQQhhlIqCEEEJoiwgoIYQQ2iICSgghhLaIgBJCCKEtIqCEUUfSNZJOGOn7CCEMTQSUEEIIbREBJYwYSdPyPiKX5/03bpW0cV2aQ/Pii4slXSVpUj7/mKRPSlqU39t9ZJ4ihFATASWMtF2AS22/GfgDcHztjbyx0zXASbbfSlrZ4R8LeVfY3g+4DPhY1+44hNBQBJQw0h61XVtmZiEwrfDebvn9h/LruaQ9SGpuapIvhDACIqCEkbay8HM/664v12ibgkZ56/OFEEZABJQwmi0Fpkl6U359KnDHCN5PCKGFCChh1LL9KmkF3/+UtJi04+NXR/auQgjNxGrDIYQQ2iJqKCGEENoiAkoIIYS2iIASQgihLSKghBBCaIsIKCGEENoiAkoIIYS2iIASQgihLf4v+t6PCGiVmNcAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sigma_z_clim = (dsz_clim.PD[k:k+2, :, :].mean('z_t') - 1.) * 1000.\n",
"sigma_z_clim = sigma_z_clim.compute()\n",
"sigma_z_clim.plot() #vmin=26., vmax=28.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the time-varing density anomaly at `z = h`"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2ae0666542e8>"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEGCAYAAABCa2PoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7xcZX3v8c83OyRAuCo3CbQJchMVCiJeOLUFVODgIa1ijVaLypFTy83bsWArrXiol1otVrRGLlJEkQZooy8EVBSrIiQgguGiaUAIFyWCgAJJ9t7f88daO0yGmb3X3pmZvfes7/v1Wq/MPPOs9axZgfnlucs2ERERG2vGZN9ARET0hwSUiIjoiASUiIjoiASUiIjoiASUiIjoiJmTfQONZm46x7O2fFb3C1L3ixgxNLtHBQ30cLReD59fL83cZKhnZW0xc01Pytlq4KmelAOw9YzhnpV14y1rVtvefmOucfghc/zrh6v9nd94y5qrbB+xMeXVwZQKKLO2fBZ7v/Y9XS9neFbXi1jv8fm9+aEf3HqwJ+UAaGYvg1fvytphh8d6VtbBO63sSTmHbnVbT8oBOGrz3gWvgef8/Bcbe41fPzzEDVf9XtXyttvY8upgSgWUiIheMTBM72pVdZCAEhG1ZMw6966Zsw4SUCKitlJD6awElIioJWOGsvRURyWgRERtDZOA0kkJKBFRSwaGElA6KgElImorNZTOSkCJiFoysC59KB2VpVciopaMGap4VCHpCEl3Sloh6dQWn8+W9NXy8+slzWv47LQy/U5JhzekbyNpsaQ7JN0u6WVN13yfJEuaEhMvU0OJiHoyDHWogiJpADgbeBWwClgqaYntxqUKjgMesb27pIXAx4A3SNoHWAg8H9gZ+JakPW0PAWcBV9o+RtIsYPOGMncty7unM99i46WGEhG1VMyUr3ZUcBCwwvZK22uBi4EFTXkWABeUrxcDh0lSmX6x7TW27wJWAAdJ2gp4BXAugO21tn/TcL1PAe8vv8qUkIASETUlhioewHaSljUcxzddbC5wb8P7VWVayzy2B4FHgWePcu5uwEPA+ZJ+LOkcSXMAJB0N3Gf7J514Ep2SJq+IqKWiU77y0tmrbR84yuetLtRcc2iXp136TOAA4CTb10s6CzhV0keAvwFePfZt91ZqKBFRS8U8lMo1lLGsAnZteL8LcH+7PJJmAlsDD49y7ipgle3ry/TFFAHmucB84CeS7i7z3yRpp6rfvVsSUCKitoatSkcFS4E9JM0vO88XAkua8iwBji1fHwNcY9tl+sJyFNh8YA/gBtsPAvdK2qs85zDgNtu32t7B9jzb8ygCzwFl/kmVJq+IqKWRGkpHrmUPSjoRuAoYAM6zvVzSGcAy20soOtcvlLSComaysDx3uaRLgNuAQeCEcoQXwEnARWWQWgm8rSM33CUJKBFRS0YMdbCRxvYVwBVNaac3vH4KeH2bc88EzmyRfjMwWt8NZS1lSkhAiYjaqticFRUloERELRmx1gOTfRt9JQElImqpmNiYcUmdlIASEbXVqU75KEytgCIYntX9YgY37X4Z602ZRRGmqR62cf/6N1v0rKyVW/RmLb+9N9u2J+UAPDj0856V1Qm2GHJqKJ00tQJKREQPDaeG0lEJKBFRS0WnfH4COylPMyJqKZ3ynZeAEhG1NZR5KB2VgBIRtdTpmfKRgBIRNTacUV4d1dWnKendkpZL+qmkr0jq5YDdiIi2isUhZ1Q6opquPSlJc4GTgQNtv4BiBc6F3SovImI8jFjngUpHVNPtJq+ZwGaS1gGb88wNZyIiJoVNJjZ2WNeepu37gE8A9wAPAI/avro5n6TjR/ZpHnzyd926nYiIJmK44hHVdLPJa1tgAcVWlTsDcyS9uTmf7UW2D7R94MzN5nTrdiIiNmCKGkqVI6rp5pN6JXCX7YdsrwMuA17exfIiIsYlnfKd1c0+lHuAl0raHHiSYj/kZV0sLyKiMlN5v/ioqGsBxfb1khYDN1Hsk/xjYFG3youIGA8D67KWFwCSDqiQbZ3tW0fL0NWnafvvgL/rZhkREROj7IfytGuBpTDqA5kPzBvtIgnPEVFLJjPlGyy1fehoGSRdM9ZFElAiorZSQymMFUyq5klAiYhaspUaSgNJOwHYflDS9sAfAnfaXl71GnmaEVFLRad8ll4BkPR/gOuAH0l6J/B14DXAZZKOq3qd1FAioqayp3yDE4HnA5sBvwB2L2sq2wLfAc6tcpGpF1DcgyJ6+K03+W1v2mg1tElPygEYeKpnRdHL/9/XbT3cs7JuG9ipJ+XMmjHUk3IANp+xpmdlFas5bZyiUz59KKV1tp8AnpD037YfBLD9iKTKv8pTL6BERPRIZsGvNyxpk3JVk6NGEsstRyo/pASUiKilzJTfwGsp24dsr2pIfzbw3qoXSUCJiNoaTg0FANv3tEm/D7iv6nXyNCOilmxYNzyj0lGFpCMk3SlphaRTW3w+W9JXy8+vlzSv4bPTyvQ7JR3ekL6NpMWS7pB0u6SXlen/WKbdIulySdts9AMprvvWjTk/ASUiaqlo8ppR6RiLpAHgbOBIYB/gjZL2acp2HPCI7d2BTwEfK8/dh2I32+cDRwCfLa8HcBZwpe29gf2A28v0bwIvsL0v8DPgtAk/iKe/w+kUq8RPWAJKRNTWULme11hHBQcBK2yvtL0WuJhiP6hGC4ALyteLgcMkqUy/2PYa23cBK4CDJG0FvIJyyK7ttbZ/U76+2vZgea0fAbtM+CEAkhYBzwP+YmOukz6UiKilcQ4b3k5S4/Ybi2w3rp4+F7i34f0q4CVN11ifx/agpEcpOr3nUgSFxnPnUmz78RBwvqT9gBuBU2w3b237duCrVb9IG28CDrK9UWPnE1AioqbGtfTKatsHjnqxZ2qev9EuT7v0mcABwEnldiBnAacCH1x/QelvKLYHuWiUe0PS0RS1HYBrbX+tKctrgEskLbD936NdazRp8oqI2urgnvKrgF0b3u8C3N8uj6SZwNbAw6OcuwpYZfv6Mn0xRYChvMaxFIHgz223nXwo6SPAKcBt5XFymbae7e9S9ON8aeyv2l4CSkTUUjHKa6DSUcFSYA9J8yXNovhxXtKUZwlwbPn6GOCaMhAsARaWo8DmA3sAN5Sz1e+VtFd5zmEUAQFJRwB/DRxdznAfzVHAq2yfZ/s8io7/o5oz2f4p8LoqX7adNHlFRC11cmJj2SdyInAVMACcZ3u5pDOAZbaXUHSuXyhpBUXNZGF57nJJl1AEi0HgBNsja+acBFxUBqmVwNvK9M8As4FvFv36/Mj2X45yi9uUZUJRM3qGcmTZGyhGoE1IAkpE1FbF5qxKbF8BXNGUdnrD66eA17c590zgzBbpNwPP6Lsphx5X9RHgx5K+Q9Ff8wrgAy2uOSRpAQkoERHjU5fFIW1/RdJ3gRdTBJS/Hln8sYUfSPoMxaix9aPJbN9UpawElIiorTpssCXp27YPo6FPpyGt2cvLP89oSDMw5m6NkIASETVli8E+DijlSsGbU8yh2ZanhydvBezc6hzbh2xMmQkoEVFbfd7k9X+Ad1EEjxt5OqA8RrFMzDNI2hH4B2Bn20eWy8K8zHalDbb6NzxHRIxipA+lyjEd2T7L9nzgfbZ3sz2/PPaz/Zk2p32RYqTaSA3mZxRBqZIElIiorX4OKCNs/8s4sm9n+xJguDx3EKi87WeavCKilrLBVku/k/RsymVjJL0UeLTqyQkoEVFbnZyHMtVIOtj2DyTNtr2m4mnvoRgN9lxJPwC2p83cmVYSUCKilmwYrLh51jT1aeBFwHU0rAE2huXAHwF7UXTi38l03VNew7DJE23XOOuYgbW9+1fJ8Ca9Kaenox+7/1e03ro5vStreHbv/rtY++jsnpTzwLZb9qQcgFuf2HXsTB1z/dhZKujzJq91ks4H5kr6dPOHtk9ucc51tg+gCCwASLqJigFpSgWUiIheqUEfymsodmA8lGLYcFuSdqLYg2UzSfuz4ZyVzasWmIASEbXlPg4otlcDF0u63fZPxsh+OPBWiqXzP9mQ/hgt1v1qJwElImqrnzvlGzwp6dvAjrZfIGlfimXv/99IBtsXABdIep3tSydaUF/3SEVEtGPXYx4K8AXgNGAdgO1bKJfOb+EHks6V9A0ASftIOq5qQQkoEVFTYmh4RqVjmtvc9g1NaYNt8p5PZspHRIyfrUrHNLda0nN5erLiMcADbfJO3ZnykrYBzgFeQPFl3m77um6WGRFRRV32QwFOABYBe0u6D7gLeHObvFN6pvxZwJW2jym3sKw8/Cwioqtc9KP0O9srgVdKmgPMsP34KNlbzZQ/pmpZXQsokrai2GryrQC21wJru1VeRMR49fsor3Kf+G1tr7b9O0mzJL0DeI/t5zXnt32TpA1mytteV7W8btZQdgMeAs6XtB/FxJpTbP+uMZOk44HjAWbN2baLtxMR8TSXnfL9StJC4PMUzVg/B/4euBBYCvx5m3MGgP8JzKOID6+WhO1PtsrfrJtPcybFdP3P2d6fYn/iU5sz2V5k+0DbB87ctIfrbERE7dnVjmnqb4EX2d4ZeDdwJXCS7T8dZY/4r1G0Kj0b2LLhqKSbNZRVwCrbI4vuLKZFQImImCx9MIJrNGttr4D1TVl32b58jHN2sb3vRAvsWkCx/aCkeyXtZftO4DDgtm6VFxExHkXto68Dyg6S3tPwfovG922asb4h6dW2r55Igd0e5XUScFE5wmsl8LYulxcRUVmfDxv+Ahs2VzW/b+VHwOWSZlDMrBdg21tVKbCrAcX2zcCB3SwjImKipnH/yJhsf2gCp/0T8DLgVnv8TyeLQ0ZELRkx3MejvCbo58BPJxJMIAElImqsjysoE/UA8N1yccj12wZXHTacgBIR9dT/nfJAMbfEdtX1uO4qj1nlMS4JKBFRX/WooqyQtBg43/aoI20n2O+yXgJKRNRWHWoowL4U+5+cU47eOg+42PZjIxkk/bPtd0n6Gi3CrO2jqxSUgBIRtWRgeLhzAUXSERQL4g4A59j+aNPns4F/A14E/Bp4g+27y89OA46jWCr+ZNtXlektV2yX9CzgqxRLpNwN/JntR1p+z2IxyC8AX5D0CuArwKfKWsuHy8mPF5bZP7Exz2BKBZQZg2bThysvvT9hw5v0bmTHmq168y+gdVv07l9aQ+NuWZ24dZUXfdh4HuhdWb3y4MOVpg90xE82mduzsjrCQIdqKOUaWGcDr6JYJWSppCVNTUzHAY/Y3r1cZ+tjwBsk7UNRg3g+xcZW35K0Z9nv0W7F9lOBb9v+qKRTy/d/Pcq9HUUxD3AexdDgi4A/BK4A9rR9Y5n9D2yf1XT+KcC1VZ5DxsxFRG11cC2vg4AVtleWK6tfDCxoyrMAuKB8vRg4TJLK9Ittr7F9F7ACOKhhxfZzi3v1Wtu/aXGtC4A/GeXefl7m/0fb+9v+pO1f2l5Msb5Xo2NbnP/W0b54oylVQ4mI6KnqnfLbSVrW8H6R7UUN7+cC9za8XwW8pOka6/PYHpT0KMUijHMpZqg3njsXeJL2K7bvaPuB8loPSNphlHv/C9vfb0yQdLDtH9g+uXz/RuBNwHxJSxqybknRPFdJAkpE1NS4tvddbXu0VT9aXag5XLXL0y59ZMX2k2xfL+ksiqatD1a430afLq/T6F+a0n5IMQdlO4omsRGPA7dULSgBJSLqq3PDhlcBuza83wW4v02eVZJmAlsDD49y7mgrtv9S0nPK2slzgF8135CklwEvB7ZvWiRyK4qBA+vZ/gXwC4plVyYsfSgRUU8GD6vSUcFSYA9J88vO84UUW+k2WsLTfRTHANeUS5wsARZKmi1pPrAHcIPtB4F7Je1VntO4YnvjtY4F/rPFPc0CtqCoODTubfIYbbb1lfRaST+X9KikxyQ9LumxVnlbSQ0lImqsM6O8yj6RE4GrKP71f57t5ZLOAJbZXkLRuX6hpBUUNZOF5bnLJV1CESwGgRMaZra3W7H9o8Alko4D7gFe3+KergWulfTFsgZSxceB/2X79vE+A0hAiYg66+BMedtXUAzDbUw7veH1U7T44S8/OxM4s0V6yxXbbf+aosbS1shkReAzkqpOVvzlRIMJJKBERJ3199IrE5msuEzSV4H/YMPFIS+rcnICSkTUUwcnNk5FI5MVy6YvACRtC+xqu93Ira2AJ4BXN14KSECJiBhNP2+wNULSd4GjKX7vbwYeknSt7fc057W9UbvqZpRXRNTXsKod09vW5UKQr6VYcfhFwCtbZZS0p6RvS/pp+X5fSX9btaAElIioLbnaMc3NLOeq/Bnw9THyfgE4jWI/ecqmsYVVC0pAiYh68jiO6e0MiuHMK2wvlbQbxfperWxu+4amtMGqBVXqQ5E02/aasdIiIqYP9XWn/Ajb/w78e8P7lcDr2mRfLem5lGFU0jEUS7JUUrVT/jqeuRZMq7SIiOlj+tc+xiRpe+AdFEvXr//Nt/32FtlPABYBe0u6j2I74DdXLWvUgCJpJ4pVLzeTtD9PTyvdiqfX5Y+ImJ6GJ/sGeuI/gf8CvkWxgVdbZe3llZLmADPKzbkqG6uGcjjFWvi7AJ9sSH8c+MB4CoqImFL6fB5Kg81tt9x8q5mkfwA+PrLvSjlv5b22K430GrVT3vYFtg8B3mr7kIbj6KozJyMipqqajPL6uqT/WTHvkQ2beFFuK1z13Gp9KLYvlXQUxRaVmzakn1G1oIiIKWf6B4sqTgE+IGktsJai68K2W+0PPdA44ErSZsDsqgVVHeX1rxR9JocA51Asfdw8tCwiIqYY21uOI/uXgG9LOp8i3L6dp7caHlPVUV4vt72vpFtsf0jSP1FxbZeIiKmqD5qzxlTuW//nwHzbH5a0K/CcFvNNsP1xSbdSrGQs4MO2r6paVtWA8mT55xOSdqbYY3h+1UKq8gwxuFn351qu3bJ38znXbdGbctZu3ZtyANbN6d3/hcOze1eWhnrXQTvz0YGxM3XA0FOb9aQcgJ+t3bFnZXWE6YdlVar4LMV4tkOBDwO/Bc4GXjySQdJVwJXAN2x/A/jGRAqq+sv6dUnbAP8I3ATcDVw8kQIjIqaMesyUf4ntE4CnYH1H+6ymPMcCjwB/L+kmSZ+TtEDSuP5JXLVT/sPly0slfR3Y1Paj4ykoImKqqUOTF7BO0gBPz37fnqYZOOV2w18EvihpBvAS4Ejg/ZKeBK62/fGxChprYuNrR/ms8qYrERFTUj0CyqeBy4EdJJ1JMaiq7bwS28MUK6FcB5wuaTuKOYljGquG8r+ayyr/FOPYdCUiYkqqQUCxfZGkG3m6o/1P2m3zK2k+cDLw+2y4TEur7YKfYdSAMrLZiqRNKRYTm9dwTg3+KiKiX/XJpMUxSXohsDfwK+D2MfaM/w/gXGAJE1iYpuoor/8AfkPRIf9UmVaDv4qI6Gt9PMpL0tYU63jtCtxCUTt5oaR7gAXlplvNnrL96YmWWTWg7GL7iIkUUHYGLQPus/2aiVwjIqIb+ryG8mGK395Dy36Rkd/jjwBnAie1OOcsSX8HXA2s357E9k1VCqwaUH4o6YW2b62Yv9EpwO0UKxRHREwd/R1QXgnsOxJMAGwPSfoA0O63/IXAWyjmrIyc5/L9mKoGlP8BvFXSXRRRa2QtmH1HO0nSLsBRFNHwPRXLiojovv7vQ1lr+xm7LdoelNRuc8Q/BXazvXYiBVYNKEdO5OLAPwPvB9quJSPpeOB4gFmbbzvBYiIiJqC/A8qmTftYjRDtF3z8CbANRQf+uFWd2PiL8V5Y0muAX9m+UdIfj3LtRRQ7hLHFs3bt77/eiJhS1N8bbD3AhvtYNXqwTfqOwB2SlrJhH8rGDxveSAcDR5fr8G8KbCXpS7YrbycZERETU+5lNV5/tzFldi2g2D4NOA2grKG8L8EkIqaUtIlswPa1G3N+75bdjYiYSiru1tjnHfcbkPRSSUsl/VbSWklDklrNV2mpJwHF9nczByUippwOrjYs6QhJd0paIenUFp/PlvTV8vPrJc1r+Oy0Mv1OSYc3pN8t6VZJN0ta1pD+B5J+NJIu6aAJfPtWPgO8Efg5sBnwv8u0SrrZhxIRMbV1qPZRThg8G3gVsApYKmmJ7dsash0HPGJ7d0kLgY8Bb5C0D7CQYov1nYFvSdrT9lB53iG2VzcV+XHgQ7a/UfZTfxz446Z7OmC0e243WdH2CkkDZfnnS/rhmA+glIASEbUkOjrK6yBghe2VAJIuBhYAjQFlAfD35evFwGfK3RQXABeX+7jfJWlFeb3rRinPPD1ZfGvg/hZ5/mmM81tNVnxC0izgZkkfpxgpNmeU62wgASUi6qmz/SNzgXsb3q+i2FOkZZ5ycuGjwLPL9B81nTv36bvkakkGPl9OswB4F3CVpE9QdF28vPmGJjjK6y3l9U4E3k2xDtjrqp6cgBIR9VU9oGzX2IcBLGr4cYdnTh5sdfV2eUY792Db90vaAfimpDtsfw94J/Bu25dK+jOKFYJf2e7mJb0A2IdiCkdRgP1vTXkGgDPL0bhPAR9qd712ElAior6qB5TVtg8c5fNVFP+aH7ELz2yGGsmzStJMiqaqh0c71/bIn7+SdDlFU9j3KLbsPaXM/+/AOe1urFzs8Y8pAsoVFCuffB/YIKCU63xtL2lWt5de6QnPgMFNu7+c9Nq2C8F03nDzzs3dKmegN+UAqKdLfvduzOaMdT0rioGnevMM9dve/V0NPtluNY+pq4NNXkuBPcoNqu6j6GR/U1OeJRSB4DqKXROvsW1JS4AvS/okRaf8HsANkuYAM2w/Xr5+NXBGea37gT8CvkvRF/LzUe7tGGA/4Me23yZpR9oHoLuBH5T39LuRRNvtZtxvYEoFlIiInupQQCn7RE4ErgIGgPNsL5d0BrDM9hKKZqkLy073hymCDmW+Syg68AeBE8rawo7A5UW/PTOBL9u+sizyHRRLzc+kaJ46fpTbe9L2sKRBSVtRrNO1W5u895fHDEZZg7GdBJSIqCd3di0v21dQNCk1pp3e8Pop4PVtzj2TYlX2xrSVFDWLVvm/D7yo4q0tk7QN8AXgRuC3wA1trjvufpNGCSgRUV81mAVv+6/Kl/8q6UpgK9u3tMor6Ws886k8SrFR1+fLoNhWll6JiNrq56VXJO1d/nnAyAE8C5g5yqTHlRQ1mC+Ux2PAL4E9y/ejSg0lIuprmgaLit5D0bfSaoJju4mN+9t+RcP7r0n6nu1XSFo+VoEJKBFRT+NYp2s6sj3SUX9kc1OVpE1bnAKwvaTfs31Pme/3gO3Kz8YcSpyAEhG1JKZvc9Y4/RBobuJqlQbwXuD7kv6b4hHNB/6qHLZ8wVgFJaBERG31c0CRtBPFEi6bNW0FvBWweatzbF8haQ9g7zL/HQ21m38eq8wElIiorz4OKMDhwFspZt43Tkx8HPhAY0ZJB4ysPlwuUvmT5os15mknASUi6quPA4rtC4ALJL3O9qVjZD+/3Fl3tKUVzgX2H+0iCSgRUU/TeEjwOH1d0puAeTT85ts+oyHP1hSTHkcLKA+NVVACSkTUVz0Cyn9STE68EVjTKoPteZ0oKAElImqrk0uvTGG72D6iFwVlpnxE1FY/z5Rv8ENJL+xFQamhREQ99fnExgb/A3irpLsomrwE2Pa+nS4oASUi6qseAeXIqhnLPe7/HNjN9hnlTPmdbLdcnbhZmrwiopZGZsr3e5OX7V9Q7Ah5aPn6Cdr/9n8WeBnwxvL948DZVctKDSUiakvD0zxaVFBuAXwgsBdwPrAJ8CXg4BbZX2L7AEk/BrD9iKTK+86mhhIR9eRxHNPbnwJHU27pW+5T3243xnWSBii/taTtgcpj4RJQIqK26tDkBay1vT40lgs9tvNp4HJgB0lnAt8H/qFqQWnyioj6mv7BoopLJH0e2EbSO4C3A+e0ymj7Ikk3AodRdDP9ie3bqxY0tQKKwAPdL2bGUPfLGDG4SW/K6eV3mvFE78qa9VjvKtEzWs4h7o6BMXeW6IwZg70pB2DNNqOt2jE19UHtY0y2PyHpVRS7L+4FnG77m63ySnopsNz22eX7LSW9xPb1VcpKk1dE1FcN+lAkfcz2N23/X9vvs/1NSR9rk/1zFFsAj/hdmVZJAkpE1JOLpVeqHNPcq1qktZuborK/BQDbw4yjJSsBJSJqqd/noUh6p6Rbgb0k3VIet5Yz5m9pc9pKSSdL2qQ8TgFWVi1zavWhRET0kqdptKjmy8A3gI8ApzakP2774Tbn/CXFSK+/pWjs+zZwfJu8z5CAEhG1NV1rH1XYfpRi2fo3StoP+MPyo/8CWgYU278CFk60zASUiKinPuhwr0LSyRS1jMvKpC9JWmT7X1rk3RQ4Dng+sOlIuu23Vymra30oknaV9B1Jt0taXrbFRURMGTXplP/fFEuqnG77dOClwDva5L0Q2IliP/prKfajf7xqQd3slB8E3mv7eRRf4ARJ+3SxvIiIcalJQBHQOFNtiPZb/e5u+4PA78o96Y8CKu+l0rUmL9sPAA+Urx+XdDswF7itW2VGRFRm+r1TfsT5wPWSLi/f/wlwbpu868o/fyPpBcCDFHvRV9KTPhRJ84D9gWfMtpR0POUogllztu3F7UREAP3dKT/C9iclfZdioy0Bb7P94zbZF0nalmKU1xJgC+CDVcvq+jwUSVsAlwLvsv1Y8+e2F9k+0PaBMzcbbc2yiIgO6+BMeUlHSLpT0gpJp7b4fLakr5afX1/+Q3vks9PK9DslHd6Qfnc5d+RmScuarndSmX+5pI+3KG9TSe+S9BngxcBnbZ/VKpg09HHfbvsR29+zvZvtHWx/vtoT6HINRdImFMHkItuXjZU/IqJXRiY2duRaxZLvZ1PMSl8FLJW0xHZjE/9xwCO2d5e0EPgY8Iayb3khxciqnYFvSdrT9ki/xyG2VzeVdwiwANjX9hpJO7S4rQsomrD+i2Jm/POAd7X5Cm8DzgL+BThgnF9/va4FlHIryXMpIt4nu1VORMSE2J3cYOsgYIXtlQCSLqb4wW8MKAuAvy9fLwY+U/5OLgAutr0GuEvSivJ6141S3juBj5bnjMwfabaP7ReW93MuMNo2vrdLupti2frGWfTj2n++mzWUgy8j57MAAArESURBVIG3ALdKurlM+4DtK7pYZkREddXjyXZNTU6LbC9qeD8XuLfh/SrgJU3XWJ/H9qCkR4Fnl+k/ajp3bsMdXi3JwOcbytwT+MNyz5KngPfZXtpU3kgH+0h5bb+c7TdK2gm4imIzrgnp5iiv79N+aFpExKQbR5PXatsHjnapFmnNV2+XZ7RzD7Z9f9mk9U1Jd9j+HsVv97YUUzJeTLHnyW6NCzsC+0ka6bcWsFn5fqTWsVVTmQ8Bt5b7zk9IZspHRD0Z6FyT1ypg14b3uwD3t8mzStJMYGuKJVDanltu14vtX5XDfg8Cvleec1kZQG6QNAxsRxEUKM8Z1+5StockbSdplu0J7diT1YYjor46N8prKbCHpPmSZlF0si9pyrMEOLZ8fQxwTRkQlgALy1Fg84E9KILEHElbwvpte18N/LQ8/z+AQ8vP9gRmARt03E/QL4AfSPqgpPeMHFVPTg0lImqrU6O8yj6KEyn6IAaA82wvl3QGsMz2EopBSheWne4PUy7CWOa7hKIDfxA4oawt7AhcXvZ9zAS+bPvKssjzgPMk/RRYCxzb1Nw1UfeXxwxgy/GenIASEbXVwVFelAOOrmhKO73h9VPA69uceyZwZlPaSmC/NvnXAm/eyFtudd0Pbcz5CSgRUU81WW14PCR9hxZPxfahVc5PQImIWiomNiaiNHlfw+tNgddRNMNVMrUCyjAMTGhswfgMrRs7T6fM/G2PylnTm3IAVPk/r403c03v/off9OGhsTN1iIZ687080LuR+3MemIazBKb/SsIdZfvGpqQfSLq26vlTK6BERPRQaigbkvSshrczgBdR7I9SSQJKRNRT+lBauZGnJ1sOAndRrEFWSQJKRNRUR9fy6gu252/M+ZnYGBH1ZVc7+pykF5dreY28/wtJ/ynp003NYKNKQImIenJttgCu4vMUEySR9Argo8C/AY8Ci0Y5bwNp8oqI+qpB7aOiAdsPl6/fQLGa8qXApQ2rxY8pNZSIqK8O7tg4zQ2UC1YCHAZc0/BZ5YpHaigRUVsarkd7VgVfAa6VtBp4kmKXRyTtTtHsVUkCSkTUk8nExpLtMyV9G3gOcHXDQpMzgJOqXicBJSJqSTgTGxvY/lGLtJ+N5xoJKBFRXwkoHZWAEhH1lYDSUQkoEVFP6UPpuASUiKitjPLqrASUiKipeiyr0ksJKBFRTyYBpcMSUCKivtLi1VEJKBFRW5mH0lkJKBFRXwkoHZWAEhH1ZMNQ2rw6KQElIuorNZSOmlIBZWCt2eLeNV0vxzPU9TJG9GqL0V5uZTrjycGelTXw6BM9K4snnuxZUUO/fKg3Bbl3/wL30FDPyuqYBJSOmlIBJSKiZwxkT/mOSkCJiJpyT2twdZCAEhH1ZNIp32EJKBFRX+lD6agElIiorwSUjpox2TcQETE5ysUhqxwVSDpC0p2SVkg6tcXnsyV9tfz8eknzGj47rUy/U9LhDel3S7pV0s2SlrW45vskWdJ2E3oEHdbVgDLWA46ImDQGhoerHWOQNACcDRwJ7AO8UdI+TdmOAx6xvTvwKeBj5bn7AAuB5wNHAJ8trzfiENt/YPvApjJ3BV4F3DOBb98VXQsoFR9wRMTk6VwN5SBghe2VttcCFwMLmvIsAC4oXy8GDpOkMv1i22ts3wWsKK83lk8B76cIjVNCN2soVR5wRMQkKZdeqXLAdpKWNRzHN11sLnBvw/tVZVrLPLYHgUeBZ49xroGrJd3YWKako4H7bP9kox5Bh3WzU77VQ3pJc6byIR0PMHv2Nl28nYiIBgZXn4eyurnJqUmr5Teaaw7t8ox27sG275e0A/BNSXcAy4C/AV49xj33XDdrKFUeMLYX2T7Q9oGzNpnTxduJiGgy7GrH2FYBuza83wW4v10eSTOBrYGHRzvX9sifvwIup2j5eS4wH/iJpLvL/DdJ2mkc37wruhlQqjzgiIjJ07k+lKXAHpLmS5pF0cm+pCnPEuDY8vUxwDW2XaYvLEeBzQf2AG6QNEfSlgCS5lDUSH5q+1bbO9ieZ3sexW/tAbYf3LiHsfG62eS1/gED91E84Dd1sbyIiOrsSiO4ql3Kg5JOBK4CBoDzbC+XdAawzPYS4FzgQkkrKGomC8tzl0u6BLgNGAROsD0kaUfg8qLfnpnAl21f2ZEb7pKuBZR2D7hb5UVEjFsHJzbavgK4oint9IbXTwGvb3PumcCZTWkrgf0qlDtvArfbFV2dKd/qAUdETA2enkvuT2FZeiUi6inL13dcAkpE1FeWr++oBJSIqCUDTg2loxJQIqKenA22Oi0BJSJqK53ynSVPof0AJD0E/GKcp20HrO7C7Uy2fvxe/fidoD+/11T/Tr9ve/uNuYCkKym+ZxWrbR+xMeXVwZQKKBMhadkYa+xMS/34vfrxO0F/fq9+/E7RfdlgKyIiOiIBJSIiOqIfAsqiyb6BLunH79WP3wn683v143eKLpv2fSgRETE19EMNJSIipoAElIiI6IhpG1AkHSHpTkkrJJ062ffTCZJ2lfQdSbdLWi7plMm+p06RNCDpx5K+Ptn30imStpG0WNId5d/Zyyb7njpB0rvL//5+Kukrkjad7HuK6WFaBhRJA8DZwJHAPsAbJe0zuXfVEYPAe20/D3gpcEKffC+AU4DbJ/smOuws4Erbe1PsWzHtv5+kucDJwIG2X0Cxl9HCyb2rmC6mZUCh2Fd5he2VttcCFwMLJvmeNprtB2zfVL5+nOIHau7k3tXGk7QLcBRwzmTfS6dI2gp4BcUufNhea/s3k3tXHTMT2Kzc93xzsnV3VDRdA8pc4N6G96vogx/eRpLmAfsD10/unXTEPwPvB/ppJb7dgIeA88umvHPKfb+nNdv3AZ8A7gEeAB61ffXk3lVMF9M1oKhFWt+Mf5a0BXAp8C7bj032/WwMSa8BfmX7xsm+lw6bCRwAfM72/sDvgGnflydpW4ra/nxgZ2COpDdP7l3FdDFdA8oqYNeG97vQJ9VySZtQBJOLbF822ffTAQcDR0u6m6Jp8lBJX5rcW+qIVcAq2yM1yMUUAWa6eyVwl+2HbK8DLgNePsn3FNPEdA0oS4E9JM2XNIui03DJJN/TRpMkijb5221/crLvpxNsn2Z7F9vzKP6errE97f/Fa/tB4F5Je5VJhwG3TeItdco9wEslbV7+93gYfTDYIHpjWu6HYntQ0onAVRSjUM6zvXySb6sTDgbeAtwq6eYy7QO2r5jEe4r2TgIuKv9RsxJ42yTfz0azfb2kxcBNFKMOf0yWYYmKsvRKRER0xHRt8oqIiCkmASUiIjoiASUiIjoiASUiIjoiASUiIjoiASWmHElflHTMZN9HRIxPAkpERHREAkpMGknzyn1EvlDuv3G1pM2a8hxWLr54q6TzJM0u0++W9CFJN5Wf7T053yIiRiSgxGTbAzjb9vOB3wCvG/mg3Njpi8AbbL+QYmWHdzacu9r2AcDngPf17I4joqUElJhsd9keWWbmRmBew2d7lZ//rHx/AcUeJCMua3NeREyCBJSYbGsaXg+x4fpyrbYpaHVu83kRMQkSUGIquwOYJ2n38v1bgGsn8X4iYhQJKDFl2X6KYgXff5d0K8WOj/86uXcVEe1kteGIiOiI1FAiIqIjElAiIqIjElAiIqIjElAiIqIjElAiIqIjElAiIqIjElAiIqIj/j/eHNEXD5jjeAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sigma_z_t = (dsz.PD[:, k:k+2, :, :].mean('z_t') - 1.) * 1000.\n",
"sigma_z_t = sigma_z_t.compute()\n",
"sigma_z_t.std('time').plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute `h_prime`, the depth of `sigma_z_t` at each time `t`"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2ae0674ca358>"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAEGCAYAAACjLLT8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZf0lEQVR4nO3df4xd5X3n8ffHHmyDgZjUwILtrUFxSMkPgusSCLtpwFl+hYXVJtmyuwmEVrJaeSmJWjXQ/FGlu0irJkpLpJaulyQlKhRRwFsUEQKbNl1Vu0Bs4vLLTrAMwRObGpNCAIN/zHz2j3sm3Fgzd8547jn3nLmfl3V0zz33nOd5jj3+znOf8/yQbSIiotnmDboAERExvQTriIgWSLCOiGiBBOuIiBZIsI6IaIGRQReg2/zjFnvkxCWDLkZfHTUyVks+I/PGa8kHwKi+vGrsrLRw/qHa8lo8b38t+Rw972At+QAsVn0/F5sf37/X9omzSePiCxb7pZ+U+/+5+fH937Z9yWzym61GBeuRE5ew7L+tH3Qx+urUk16uJZ8TFr5RSz4A467vP+WbY/X9iK46/sXa8jrnuGdryee9i0ZryQfglxcsqC2v+ac886PZpvHST8Z49Nv/smx+S2eb32w1KlhHRNTFwDj1fSOdrQTriBhKxhx0Pc2U/ZAHjBExtMZL/ilD0hJJd0vaJmmrpPMkfbF4/7ikjZKWdJ1/o6Ttkn4g6eLp0k+wjoihZMyYy20l3Qw8YPtdwFnAVuAh4D223wf8ELgRQNKZwFXAu4FLgD+TNL9X4gnWETG0xnGpbTqSjgc+BHwVwPYB2y/bftD2RDejh4Hlxf6VwJ2299t+FtgOnNMrjwTriBhKBsZwqa2E04EXga9L+r6kWyUtPuycXwe+VewvA3Z2fTZaHJtSgnVEDK0Z1KyXStrUta07LKkRYDVwi+2zgdeBGyY+lPR54BBw+8ShSYrT87dCeoNExFAycLB8e/Re22t6fD4KjNp+pHh/N0WwlnQNcDmw1m/NST0KrOi6fjmwq1cBUrOOiKHkkk0gZZpBbL8A7JR0RnFoLfC0pEuAzwFX2N7Xdcl9wFWSFko6DVgFPNorj9SsI2I4Gcb6O53BdcDtkhYAO4Brge8BC4GH1BmO/7Dt37T9lKS7gKfpNI+st3t3+k6wjoih1BnB2Mf07C3A4U0l7+hx/k3ATWXTT7COiCElxmqclGy2EqwjYih1HjAmWEdENFqnn3WCdURE49U53e9sJVhHxFBKzToiogWMGGvRUJME64gYWmkGiYhoOCMOuOespI2SYB0RQ6kzKCbNIBERjZcHjEdqXIzvb1aRZmvPK8fVks9rixbWkg/A+Hh9P+AHD9X3NfWlfYdPP1ydVw4eXU8+bzumlnwA9o0/V1te/WCLMadmHRHReOOpWUdENFvnAWN7QmB7ShoR0Ud5wBgR0RJj6WcdEdFsGcEYEdES4y3qDVJpSSV9VtJTkp6U9FeSFlWZX0REWZ2JnOaV2pqgslJIWgb8NrDG9nuA+cBVVeUXETETRhz0/FJbE1TdDDICHC3pIHAM0yy1HhFRF5tWDYqprKS2fwx8CXge2A28YvvBw8+TtE7SJkmbxl57variREQcRoyX3JqgymaQE4ArgdOAU4HFkj55+Hm2N9heY3vN/GPrG+4bEcPNdGrWZbYmqLIUHwGetf2i7YPAvcAHK8wvImJG8oCx43ngXEnHSBKwFthaYX4REaUZMe5yWxmSlki6W9I2SVslnSfp7ZIekvRM8XpCca4kfUXSdkmPS1o9XfpVtlk/AtwNPAY8UeS1oar8IiJmwsBBj5TaSroZeMD2u4Cz6FRObwC+Y3sV8J3iPcClwKpiWwfcMl3ilfYGsf0HwB9UmUdExJFR3+azlnQ88CHg0wC2DwAHJF0JfLg47Tbgu8Dn6DzP+4ZtAw8XtfJTbO+eKo9mNMZERNTMdEYwltmApRO91opt3WHJnQ68CHxd0vcl3SppMXDyRAAuXk8qzl8G7Oy6frQ4NqUMN4+IoTWDmvVe22t6fD4CrAaus/2IpJt5q8ljMpNl7F4FSM06IoaSrZnUrKczCowWz+qg87xuNfBPkk4BKF73dJ2/ouv65UwzaDDBOiKGUucBY3+Gm9t+Adgp6Yzi0FrgaeA+4Jri2DXA3xT79wFXF71CzqUzaHDK9mpIM0hEDK2+r8F4HXC7pAXADuBaOhXiuyT9Bp3uzJ8ozr0fuAzYDuwrzu2pWcF6TMx7ufoieaRn01BfHXy9ngVLX55fzwKsAPPfqPELWX3/VLyxeLy2vB59rZ4JKHf+wpJa8gHY8fYTa8sLnpl1Cp0HjP0bSm57CzBZu/baSc41sH4m6TcrWEdE1KgpoxPLSLCOiKE0MYKxLRKsI2JoZcHciIiGs+HgeIJ1RESjdZpBEqwjIhqvX3OD1CHBOiKGUr+77lUtwToihlSaQSIiWqEp6yuWkWAdEUOp0xtk+nk/miLBOiKGUgbFRES0RJpBIiIaLr1BIiJaIr1BIiIazhaHEqwjIpovzSAREQ2XNuuIiJZIsI6IaLj0s46IaIn0s46IaDgbDmXxgSMz7xAserH633Sq86tPTatzq76FuZl3qL68xo6qL68DJ9Q3T8Shg/Wsbr67xp/1A4caFU5KaVMzSHt+rURE9NFEm3WZrQxJz0l6QtIWSZuKY++X9PDEMUnnFMcl6SuStkt6XNLq6dJv36/CiIg+cf9r1hfY3tv1/o+AL9j+lqTLivcfBi4FVhXbB4BbitcppWYdEUNrHJXaZsHA8cX+24Bdxf6VwDfc8TCwRNIpvRJKzToihpI9ozbrpRNNG4UNtjccniTwoCQD/6P4/DPAtyV9iU7l+IPFucuAnV3XjhbHdk9VgATriBhSYqx8b5C9ttdMc875tndJOgl4SNI24OPAZ23fI+k/AF8FPgKTVtd7dkdIM0hEDC1bpbZyaXlX8boH2AicA1wD3Fuc8tfFMejUpFd0Xb6ct5pIJlVpsJa0RNLdkrZJ2irpvCrzi4goa2JukH70BpG0WNJxE/vARcCTdALwrxanXQg8U+zfB1xd9Ao5F3jF9pRNIFB9M8jNwAO2Py5pAXBMxflFRJTjTrt1n5wMbJQEnbh6h+0HJL0G3CxpBHgTWFecfz9wGbAd2AdcO10GlQVrSccDHwI+DWD7AHCgqvwiImaqX8PNbe8Azprk+D8AvzzJcQPrZ5JHlTXr04EXga9LOgvYDFxv+/XukySto/htM3L8CRUWJyLiLZ7ZA8aBq7KkI8Bq4BbbZwOvAzccfpLtDbbX2F4zcsziCosTEfHz7HJbE1QZrEeBUduPFO/vphO8IyIaoZ+9QapWWbC2/QKwU9IZxaG1wNNV5RcRMROdWnN7gnXVvUGuA24veoLsoMQTz4iIurRp1r1Kg7XtLcB0o34iIgaiKe3RZWS4eUQMJSPGW9QbJME6IoZWiyrWCdYRMaRcyXzWlUmwjojh1aKqdYJ1RAyt1KwjIhrOwPh4gvURmXcQFr9Q/feSOn+Zjh9VT2bjNf5Lur5FwGFBfVmN7KsvL43V0wvh0KF6VlEH2PNGo8LJ9Ey9wWCWWva3GxHRP+lnHRHRBgnWERFN15x5P8pIsI6I4ZWadUREwxmc3iAREW2QYB0R0XxpBomIaIEE64iIhmvZoJj2TOYaEdFn/VwwV9Jzkp6QtEXSpq7j10n6gaSnJP1R1/EbJW0vPrt4uvRTs46I4dX/3iAX2N478UbSBcCVwPts75d0UnH8TOAq4N3AqcD/lvRO22NTJZyadUQMLbncNgu/Bfx32/sBbO8pjl8J3Gl7v+1nge3AOb0SSrCOiOHkGWzlU3xQ0mZJ64pj7wT+taRHJP29pF8pji8DdnZdO1ocm1KpZhBJCyd+M/Q6FhHRHprJA8al3e3QwAbbGw4753zbu4qmjockbaMTY08AzgV+BbhL0ulM3sG756+Fsm3W/w9YXeJYRER7lK8177W9pmdS9q7idY+kjXSaNUaBe20beFTSOLC0OL6i6/LlwK5e6fcM1pL+BZ2q+dGSzuat3wbHA8f0ujYiovHG+5OMpMXAPNuvFvsXAX8IvAZcCHxX0jvpzNC+F7gPuEPSl+k8YFwFPNorj+lq1hcDn6YT9b/cdfxV4PdnekMREY3R337WJwMbJUEnrt5h+wFJC4CvSXoSOABcU9Syn5J0F/A0cAhY36snyESiU7J9G3CbpI/Zvmf29xMR0Ryz7OnxM7Z3AGdNcvwA8MkprrkJuKlsHqXarG3fI+mjdPoELuo6/odlM4qIaJwWDTcv1XVP0p8DvwZcR6fd+hPAL1ZYroiI6FK2n/UHbV8N/LPtLwDn8fNPMiMiWqeGQTF9U7br3hvF6z5JpwIvAaf1uzAyzN9f/d/MoUU1Tt5S0z+0axzeVOdK6nV+TZ3/Zn15jbxWTz5H/bS+n/Wxnx5VW159YaoYbl6Zsv/tvilpCfBF4DE6t3lrZaWKiKhDQ2rNZZR9wPhfi917JH0TWGT7leqKFRFRvaY0cZQx3aCYf9/jM2zf2/8iRUTUZK4Ea+DfHvZ+4tZU7CdYR0R7zZVgbftaAEmLgI8BK7uuadFtRkT8vCb19Cij7APG/wW8TOfh4sQz8xbdZkTEJOZgb5Dlti85kgwkzQc2AT+2ffmRpBERUYU21azL9s79v5Lee4R5XA9sPcJrIyKq09/FBypVNlj/K2BzsbDj48WikI9Pd5Gk5cBHSZ/siGiakqMXm1L7LtsMcukRpv8nwO8Bx011QrH8zTqABceccITZREQcgYYE4jLKDor50UwTlnQ5sMf2Zkkf7pH2BmADwLG/sKJFf3UR0Xbq0+IDdahyRonzgSskPQfcCVwo6S8rzC8iYs6qLFjbvtH2ctsrgauAv7U96STcERED0aIHjHXOnxYR0RwNenhYRi3B2vZ3ge/WkVdERGkJ1hERLZBgHRHRbKJdvUESrCNiOKXNOiKiJVoUrGtcuS8iomH62HVP0nPFVBxbJG067LPflWRJS4v3kvQVSduLKTxWT5d+s2rWrqkNqcZZEeeN1fSr+0B9NzVvrLasGK9xDVbVeF/z99eTj+r6+QPGFrZnutEJFTSDXGB778/lIa0A/g3wfNfhS4FVxfYB4JbidUqpWUfE8KpnUMwf05kjqTulK4FvuONhYImkU3olkmAdEcOp+CZfZgOWStrUta2bPEUelLR54nNJV9CZy/8fDzt3GbCz6/1ocWxKzWoGiYioU/la817ba6Y553zbuySdBDwkaRvweeCiSc6drM2oZ2kSrCNiaPWzzdr2ruJ1j6SNwK8CpwH/KAlgOfCYpHPo1KRXdF2+HNjVK/00g0TE8OpTm7WkxZKOm9inU5v+nu2TbK8sJrQbBVbbfgG4D7i66BVyLvCK7d298kjNOiKGU39n1DsZ2FjUoEeAO2w/0OP8+4HLgO3APuDa6TJIsI6IoST61wxiewdw1jTnrOzaN7B+JnkkWEfE0Mpw84iINkiwjohogQTriIiGy6x7EREtkWAdEdF8WXwgIqIF0gwSEdF0/R0UU7kE64gYXgnWERHN1s8RjHVIsI6IoaXx9kTrBOuIGE5ps46IaIc0g0REtEGC9REyzDtY/d/e/P31rcJ8aFE9ec3fX99PXZ0DCerMq86/w7pWva93dfP2rWWSmnVERBskWEdENJwz3DwiovHSzzoioi3cnmidYB0RQys164iIpmvZoJjK+tpIWiHp7yRtlfSUpOuryisi4khovNzWBFV2jDwE/I7tXwLOBdZLOrPC/CIiZqSfwVrSc5KekLRF0qbi2BclbZP0uKSNkpZ0nX+jpO2SfiDp4unSryxY295t+7Fi/1VgK7CsqvwiImbEdB4wltnKu8D2+22vKd4/BLzH9vuAHwI3AhQV16uAdwOXAH8maX6vhGsZciRpJXA28Mgkn62TtEnSpoP7X6ujOBERQOcBY5ntSNl+0Pah4u3DwPJi/0rgTtv7bT8LbAfO6ZVW5cFa0rHAPcBnbP/08M9tb7C9xvaaoxYeW3VxIiLe4pIbLJ2oVBbbuilSe1DS5ik+/3XgW8X+MmBn12ejTNPyUGlvEElH0QnUt9u+t8q8IiJmYoaDYvZ2NW1M5XzbuySdBDwkaZvt/wMg6fN0nuPd3pX94XqWpsreIAK+Cmy1/eWq8omIOCI2Gi+3lUvOu4rXPcBGimYNSdcAlwP/2f5ZA/gosKLr8uXArl7pV9kMcj7wKeDC4unoFkmXVZhfRMTMlG8G6UnSYknHTewDFwFPSroE+Bxwhe19XZfcB1wlaaGk04BVwKO98qisGcT2PzB5VT8iohH6OILxZGBjp0GBEeAO2w9I2g4spNMsAvCw7d+0/ZSku4Cn6TSPrLc91iuDjGCMiOFkoE9rMNreAZw1yfF39LjmJuCmsnkkWEfE8GrRcPME64gYWpnIKSKiBcr29GiCBOuIGE4tm3UvwToihlJnUEx7onWjgrVsRmpYYdrz594q1gterW8ex3n7a8zrUH15jbx6oLa86lqhxPNqXHG8fYubQ0OmPy2jUcE6IqJOqVlHRDRd2qwjItqg/LwfTZBgHRHDK80gEREN5+asr1hGgnVEDK/UrCMiWqA9sTrBOiKGl8bb0w6SYB0Rw8lkUExERNMJZ1BMREQrJFhHRLRAgnVERMOlzToioh3SGyQiovGcZpCIiMYzrQrWbZwuPCKiP8ZLbiVIek7SE5K2SNpUHHu7pIckPVO8nlAcl6SvSNou6XFJq6dLP8E6IoaW7FLbDFxg+/221xTvbwC+Y3sV8J3iPcClwKpiWwfcMl3CCdYRMbzsctuRuxK4rdi/Dfh3Xce/4Y6HgSWSTumVUIJ1RAwnG8bGy20lUwQelLRZ0rri2Mm2d3ey827gpOL4MmBn17WjxbEp5QFjRAyv8rXmpRPt0IUNtjccds75tndJOgl4SNK2HulpstL0KkCjgvW8g+MsemFf5fksqrFvpWpa3ZyxsXryAfTmwdry4o03a8vKb9aX19hP/rmejFTjl2e3p8/yz5QP1nu72qGnSMq7itc9kjYC5wD/JOkU27uLZo49xemjwIquy5cDu3qln2aQiBhOBsZdbpuGpMWSjpvYBy4CngTuA64pTrsG+Jti/z7g6qJXyLnAKxPNJVNpVM06IqI+7ue3gZOBjZKgE1fvsP2ApO8Bd0n6DeB54BPF+fcDlwHbgX3AtdNlkGAdEcPJzOThYe+k7B3AWZMcfwlYO8lxA+tnkkeCdUQMrxaNYEywjojhlWAdEdF07ZrIqdLeIJIukfSDYvz7DdNfERFREwPj4+W2BqgsWEuaD/wpnTHwZwL/UdKZVeUXETFj1Q8375sqm0HOAbYXT0mRdCed8fBPV5hnRERJ7ltvkDpUGawnG/v+gcNPKsbQrwNYtOBtFRYnIqKLwS0adVllsC419r0YX78B4G2LT23G942IGA4lRic2RZXBesZj3yMiatWQ9ugyqgzW3wNWSToN+DFwFfCfKswvIqI8uzE9PcqoLFjbPiTpvwDfBuYDX7P9VFX5RUTMWGrWHbbvpzNhSUREwxjXOLXwbGUEY0QMp4kpUlsiwToihle67kVENJsBp2YdEdFw7uviA5VLsI6IodWmB4xyg7quSHoR+NEML1sK7K2gOIM2F+9rLt4TzM37avo9/aLtE2eTgKQH6NxnGXttXzKb/GarUcH6SEjaNN2qw200F+9rLt4TzM37mov31HZZ3TwiogUSrCMiWmAuBOsNgy5ARebifc3Fe4K5eV9z8Z5arfVt1hERw2Au1KwjIua8BOuIiBZobbCeiyunS1oh6e8kbZX0lKTrB12mfpE0X9L3JX1z0GXpF0lLJN0taVvxb3beoMvUD5I+W/z8PSnpryQtGnSZoqXBeg6vnH4I+B3bvwScC6yfI/cFcD2wddCF6LObgQdsvws4izlwf5KWAb8NrLH9Hjpz0V812FIFtDRY07Vyuu0DwMTK6a1me7ftx4r9V+n851822FLNnqTlwEeBWwddln6RdDzwIeCrALYP2H55sKXqmxHgaEkjwDFkOb5GaGuwnmzl9NYHtW6SVgJnA48MtiR98SfA7wHtmTVneqcDLwJfL5p3bpW0eNCFmi3bPwa+BDwP7AZesf3gYEsV0N5gXWrl9LaSdCxwD/AZ2z8ddHlmQ9LlwB7bmwddlj4bAVYDt9g+G3gdaP2zE0kn0PmWehpwKrBY0icHW6qA9gbrObtyuqSj6ATq223fO+jy9MH5wBWSnqPTXHWhpL8cbJH6YhQYtT3xzeduOsG77T4CPGv7RdsHgXuBDw64TEF7g/XPVk6XtIDOA5D7BlymWZMkOm2gW21/edDl6QfbN9pebnslnX+nv7Xd+pqa7ReAnZLOKA6tBZ4eYJH65XngXEnHFD+Pa5kDD07nglbOZz2HV04/H/gU8ISkLcWx3y8WHo7muQ64vagw7ACuHXB5Zs32I5LuBh6j0zvp+2ToeSNkuHlERAu0tRkkImKoJFhHRLRAgnVERAskWEdEtECCdURECyRYR+NI+gtJHx90OSKaJME6IqIFEqxjYCStLOaB/p/F/MkPSjr6sHPWFhMlPSHpa5IWFsefk/QFSY8Vn71rMHcRUY8E6xi0VcCf2n438DLwsYkPiknv/wL4NdvvpTPi9re6rt1rezVwC/C7tZU4YgASrGPQnrU9MbR+M7Cy67Mzis9/WLy/jc4c0hPuneK6iDknwToGbX/X/hg/P1/NZFPhTnbt4ddFzDkJ1tFk24CVkt5RvP8U8PcDLE/EwCRYR2PZfpPOTHZ/LekJOivN/PlgSxUxGJl1LyKiBVKzjohogQTriIgWSLCOiGiBBOuIiBZIsI6IaIEE64iIFkiwjohogf8Pn8W7XEQmvqkAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"nl = len(dss.time)\n",
"nj = len(dss.nlat)\n",
"ni = len(dss.nlon)\n",
"\n",
"h_prime = xr.DataArray(np.ones((nl, nj, ni)) * np.nan, \n",
" dims=('time', 'nlat', 'nlon'), coords={'time': dss.time})\n",
"\n",
"for j in range(len(dss.nlat)):\n",
" for i in range(len(dss.nlon)):\n",
" if not np.isnan(sigma_z_clim[j, i]):\n",
" h_prime[:, j, i] = dss.Z.isel(nlat=j, nlon=i).interp(sigma=sigma_z_clim[j, i])\n",
"\n",
"h_prime.std('time').plot()"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2ae067614b00>"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZQdxX328e+jfUOWxGaQMIIEs9oIULCwDSaAbRz7NQQvQGIQHBLFJCZektfGThy/dpZDHC+B9zhwiFlE4GWJIIEQlhBsfExMFAQIsUgggWU0IIQWJCSBkDTze/+outAMs/RI996Z2/f5nFNnbld1d1WPRr+pqa6uVkRgZmbVMmywG2BmZvXn4G5mVkEO7mZmFeTgbmZWQQ7uZmYV5OBuZlZBDu4NJuldkjZJGj7YbTGz9uHgXmeSlks6qbYdEc9FxISI6BzMdvVF0sGSfiJpg6Rlkn67UDZK0rx8XSHp+G7HStLfSlqb03clqVA+XNJfSXpB0kZJj0iaVCjfX9LtuWyNpO/20sZj8y/JYgpJnyrs82VJL+bruFLS6ELZDEk/z2Udkv6iUHaIpAWSXs7pPyUdUij/kqRnJb2Sr+OHkkYUyn8qaXUuf1TSKYWy35T0mKT1+fvzL5KmFsq/K2lFPvZXkv6sULabpP/Kx62X9ICkDxTKz5D0VL6mlyTNlTSxUH6tpJX53E9L+r1efgSsiiLCqY4JWA6cNNjtGEB7RwBPA18BhgMnAJuBd+fyUcCXgA8CK4Hjux3/B8BTwDRgKvAk8PlC+V8BPwH2BQQcBowpnPuZXPd4YAzw3pLtPh7YCIzP2x8FVgGHApOB+4CLCvs/Cfx1vsZfy9fyyVw2CZie2zcc+GNgUeHYXwMm5c9T8vV8pVD+XmBE/vy+3K698vaewN7582jgu8BthWMPLFzDVOAJ4LS8PSaXD8ttOxVYV6hrH2C3/HkCcB1wSeHchwKj8+eDgBeBowb7Z86pOck99zqS9E/Au4B/yz3Lr0qannuYI/I+9+We7C/yPv8maVdJ1+Ue1oOSphfOeZCkeySty720z9a52QcBewM/jIjOiPgJ8F/AWQARsTUi/j4i7gd6+utjNvD9iOiIiOeB7wPn5LZPJv1i+P2I+FUkj0fElnzsOcALEfGDiNgcEVsiYlHJds8G5kXE5sL2FRHxRES8DPxlrR3ZdOC6fI3PAPeTgh8RsT4ilkdEkIJoJ/DrtQMj4pmIWJ83BXR1K18UEdtrm8BIUuAlIlZFxAuFdnQ/91OFa6B47vz9eCoiugrtmkz6BUNErIiINX2c+4mIeL3QriD9orJ2MNi/XaqW6NZzJwWV4M3e1n3AMtJ/sneQepRPAyeRetHXAFflfccDK4Bzc9mRwBrg0F7q/gdgfS9pUS/HvAfYBKiQdw/wLz3s28Hbe+4bgPcVtmcCG/Pn43LdXyP1Gp8G/qiw75XAPwF35uu6D3hPie/xOFLv+PhC3qPA6YXt3fL3fde8/TfARaTAe2C+lt/odt71wHZSgP3zbmW/A7ySz7kaOLxb+e3Allx+FzCsUPaufO4uYBtwTrdjL8z/BgE8C0zrVr4I2JrL/7Fb2Qfzv0GQ/uL6SA8/E6/m8oeBCYP9f8SpOck998FxVaTe4AZSYHsmIv4zUu/vn4Ej8n6fAJZHxFURsT0iHgZuBj7d00kj4g8jYlIv6b29tGUJ8BLwvyWNlPQR4EOkAFrGBFJwqdkATMjj7tNIv8DeDeyX2/1/JH047zsNOAO4hPTXw78Dt0oa1U+dnyL9MvhZP+0A2CV/vT3X/xrpmq+IiAeLJ42ISbm9XwAe6Vb2/yJiYr6Wy0hDQMXyT+S6fgu4O1Jvu1b2XD73bsCf5/qLx16Ujz2S9MtuQ7fy9wITSb9g7u9Wdn9EvIP0vfw7UueiWP6H+dzHArcAr2NtwcF9cBQDw2s9bE/In/cF3pdvpq2XtB74XeCd9WpIRGwjjeV+nNS7/hPgJlLPtoxNpMBTMxHYFBFBuhaA70TEa5GGXG4gBUBy+f0RcWdEbAW+B+wKHNxPnbOBa3IdfbUDYKOkKaTe9HdI49j7AB+V9IfdTxxpiOQy4BpJe/RQvpQ0Lv4PPZRti4g787k/2UP5OmAu6RfYiG5lERGPkL4n3+7h2C0RcT1woaTDeyh/Pl/jDT2UdUYaVpsGnN+93KrJwb3+6rnM5grgZ9164BMiosf/oJIu62FGSS090WuD05jxhyJi14j4KLA/8D8l2/gEUAw2h+c8SMMJ0Pv3ZFEfZT2StA/pZuo1JdqxKiLWkq6nMyKuyX8BdfDWXzLdDSP95TK1l/IR9D123Vf5CGAP3vqLaCDnHkm6nnq3yyrGwb3+VtH7f76Buh14t6Sz8pDJSEm/IanHnm1EfD4H/57Sob1VIum9ksZIGifpT4G9gKsL5aMljcmbo/K+temO1wBfkTRV0t6knv/VuT3PAD8H/iyf42Dg9HxdANcCsySdpPQcwJdIwy2L+/ienAX8Ip+76BrgvDytcTJp+KN2DU+ny9DvSBom6Z25HY/m6/uwpCOUpm1OBH4AvFxrh6Tfq/Xi8xTJrwP35u2DJH1M0tj87/M50r2Gn+Xy0yQdmOvdPZ/7kYhYl/P+QNJkJUcDf1Q49yxJH1SajjpW0tdIs2/m5/LfVXqOQpL2Jc0Gqh27R54qOSFf10eBM0kzfawdDPagf9UScArwHOkG2p/S8w3V3yvs/1fA1YXtk4Blhe0DSWPRq4G1pP+cM+rc5r8jBbNNpHsAv96tfDlvzraopem5TKTpfety+i5vvTk7lTRcsIl0s/APup37NNIN5lfy9+bQQtmdwDe67b8EOK+X6/gK6ZfrK8BV5GmAuewE4EHSePaLwD8C43LZZ/J5N+Xv8x0UpmTmc60i3bBcnr9ftemcB5OC7cb8b/4g8NuFYy8AfpmPfZH0F8O+uWxY/t6sy3U/DXyj9v0j3ft4NJ97HekXxnGFc/81afhsc/56OW/eQN49778+fz8eI81aGvT/I07NSbUfIjMzqxAPy5iZVZCDu5lZBTm4m5lVkIO7mVkFjeh/l+YZMWZ8jN5lSsPrGba1iTeRm1RV1yj1v1O9NLOqJq6lOXzT1qbVFdu2NaUeje7vYd/66RrTvHCyacPzayJi9505x0d/c3ysXVfuB+yhRa/fHREn70x9zTakgvvoXaZw0KlfaXg9E15ozn8sgGHbmhPdN00b2ZR6ALqGNy+6j1nf1f9OdTJx/oqm1bW94/mm1DNin/2aUg/A5kN2KtYOyM9v++qvdvYca9d18j93v6vUvsP3WrrbztbXbEMquJuZNUsAXTSv89BsDu5m1paCYNvQfYfOTnNwN7O2VeWeu2fLmFlbCoLOKJfKkDRJ6ZWUSyQtlnSM0usd/1vSQqVXOR6d95WkS5Rea7lI0pGF88yWtDSn2YX8o5Re2bgsH9vnzS8HdzNrW11EqVTSxcBdEXEQaVXSxaS1lr4dETOAv8jbAB8DDshpDnApQF6e+luk1zUeDXwrL4RH3mdO4bg+Z+84uJtZWwqgkyiV+pNXEz0OuALeeD3l+lxNbXnndwC1Vy6eQn4nQUT8NzBJ0l6kdwHfExHrIr0u8h7g5Fw2MSIeiLQg2DWk9zD0ymPuZta2BtAr303SgsL25RFxeWF7f9KKolfll6k8BHyRtIz13ZK+R+pMvz/vP5X0voaajpzXV35HD/m9cnA3s7YUwLbyq+KuiYiZfZTX3nF8QUTMl3Qx6d247wC+HBE3K73c/grSst49jZfHDuT3ysMyZtaWouSQTJlhGVJPuiMi5ufteaRgP5v07lpI70c+urD/PoXjp5GGbPrKn9ZDfq8c3M2sPQV0lkz9niriRWCFpANz1onAk6QA/KGcdwKwNH++DTg7z5qZBWyIiJXA3cBH8tu5JgMfIb1wfSXpfcCz8iyZs4Fb+2qTh2XMrC2lJ1Tr6gLgOkmjSG8dO5cUgC/OL0TfQprtAultX79FegvZq3lfIr1+8S9Jb/SC9HL5dfnz+aRXR44lvaXszr4a4+BuZm1KdNZxFbyIWAh0H5e/Hziqh32D9L7cns5zJXBlD/kLgMPKtsfB3czaUrqh2sQlTpvMwd3M2lKa5+7gbmZWOV3uuZuZVYt77mZmFRSIzgrPBndwN7O25WEZM7OKCcTWGD7YzWgYB3cza0vpISYPy5iZVY5vqDbJ9jGw/uDSq7TtsM3TRja8jprXpm5vSj3DNzalGgDU1fh/ozfq6mxez2rDfvs2ra6Rm5pT19ZdmlINAK+9s3k/F9y286eIEJ3hnruZWeV0ueduZlYt6YZqdUNgda/MzKwPvqFqZlZRnZ7nbmZWLX5C1cysoroqPFumoVcm6cuSnpD0uKTrJY1pZH1mZmWlhcOGlUqtqGGtljQV+GNgZkQcBgwHzmhUfWZmAxGIbTG8VGpFjR6WGQGMlbQNGEc/b+s2M2uWCCr9EFPDriwinge+BzwHrCS93fs/uu8naY6kBZIWdG3e3KjmmJl1I7pKplbUyGGZycApwH7A3sB4SZ/rvl9EXB4RMyNi5rDx4xvVHDOztwhSz71MakWNbPVJwC8jYnVEbANuAd7fwPrMzAakyjdUGznm/hwwS9I44DXgRGBBA+szMystkF/WsSMiYr6kecDDwHbgEeDyRtVnZjYQAWzz2jI7JiK+BXyrkXWYme0YeT13M7OqCar9hKqDu5m1LffczcwqJkLuuZuZVU26odqaSwuU4eBuZm3K71BtmnHjX2fGrGUNr+fxF9/Z8DpqJgzvako9I/fqbEo9zfb61ub9iG6cMq5pdQ1/tTlBJaZsbUo9AHvu/krT6vplHc6Rbqh6zN3MrHJa9enTMqp7ZWZmfag9oVomlSFpkqR5kpZIWizpGEk3SlqY03JJC/O+oyRdJekxSY9KOr5wnqNy/jJJl0hSzp8i6R5JS/PXyX21x8HdzNpWF8NKpZIuBu6KiIOAw4HFEXF6RMyIiBnAzaQ1tgB+HyAi3gN8GPi+pFpFlwJzgANyOjnnXwjcGxEHAPfm7V45uJtZW4qAbV3DSqX+SJoIHAdckc4dWyNifaFcwGeB63PWIaQATUS8BKwHZkraC5gYEQ9ERADXAKfmY04B5ubPcwv5PXJwN7O2lIZlhpVKwG61907kNKfb6fYHVgNXSXpE0o8lFdcwPxZYFRFL8/ajwCmSRkjaDzgK2AeYCnQUjuvIeQB7RsRKgPx1j76uzzdUzaxtDeAJ1TURMbOP8hHAkcAFedHEi0nDJt/M5WfyZq8d4ErgYNJKub8CfkFaYLGnBkXZRnZvkJlZ26nzVMgOoCMi5ufteeQxcUkjgNNIvfNUd8R24Mu1bUm/AJYCLwPTCuedxpuvJ10laa+IWJmHb17qq0EeljGzNjWgYZk+RcSLwApJB+asE4En8+eTgCUR8cZwi6RxtWEbSR8GtkfEk3m4ZaOkWXmc/mzg1nzYbcDs/Hl2Ib9H7rmbWduq8/tRLwCukzQKeBY4N+efwVuHZCCNl98tqQt4HjirUHY+cDUwFrgzJ4CLgJsknUd6GdJn+mqMg7uZtaU0W6Z+a8tExELgbePyEXFOD3nLgQO75+eyBcBhPeSvJf1FUIqDu5m1Jb9mz8ysouo8LDOkOLibWVvywmFmZhXll3WYmVVMhNju4G5mVj0eljEzqxiPuZuZVZSDu5lZxXieu5lZRXmeu5lZxUTA9hIv4mhVQyq4v2PEq3x8t0UNr2fPMc17S3uz/uxr5nzdzZ2jmlbXpm2jm1bXC2MmNq2uza8353u458SNTakHYMaU55tW1/z+dynFwzJmZhXjMXczs4oKB3czs+rxDVUzs4qJ8Ji7mVkFiU7PljEzq54qj7k39NeWpEmS5klaImmxpGMaWZ+ZWVm1tWXKpFbU6J77xcBdEfHp/NLYcQ2uz8ysnEjj7lXVsOAuaSJwHHAOQERsBbY2qj4zs4Gq8myZRg7L7A+sBq6S9IikH0sa330nSXMkLZC0YNPL2xvYHDOzN0W+oVomtaJGtnoEcCRwaUQcAWwGLuy+U0RcHhEzI2LmhMm+v2tmzRNRLrWiRgb3DqAjImrLQMwjBXszsyEhQqVSK2pYcI+IF4EVkg7MWScCTzaqPjOzgUi98uoG90aPg1wAXJdnyjwLnNvg+szMSmvVaY5lNDS4R8RCYGYj6zAz21GtOp5ehu9gmllbCkRXi86EKcPB3czaVoU77g7uZtamotpryzi4m1n7qnDX3cHdzNpWlXvu1b2bYGbWhwC6ulQqldHTKriSbpS0MKflkhbmfUdKmivpsbzv1wvnOVnSU5KWSbqwkL+fpPmSlubz9vmW9SHVc5847HVOGPdMw+v5tVEvNbyOmqVb92xKPbsM29KUegA6m9gnWLalOd8/gCmjd21aXR2bJzWlnj3GbmpKPQCHjnu+aXXVRQD17bm/bRXciDi9Vijp+8CGvPkZYHREvEfSOOBJSdcDK4AfAR8mPeX/oKTbIuJJ4G+BH0bEDZIuA84DLu2tMe65m1nbqtfaMoVVcK9I542tEbG+UC7gs8D1taqB8ZJGAGNJK+a+AhwNLIuIZ/NKujcAp+TjTyAt4wIwFzi1rzY5uJtZ+4qSCXarrV6b05xuZ+pvFdxjgVURsTRvzyMtprgSeA74XkSsA6aSeu81HTlvV2B9RGzvlt+rITUsY2bWPANaN2ZNRPT1tH1tFdwLImK+pItJq+B+M5efyZu9dkg99E5gb2Ay8HNJ/wk9LjAffeT3yj13M2tf5Xvu/el1Fdw89HIacGNh/98hjc9vi4iXgP8iLdXSAexT2G8a8AKwBpiUz1XM75WDu5m1p4DoUqnU76n6XgX3JGBJRHQUDnkOOEHJeGAWsAR4EDggz4wZBZwB3BYRAfwU+HQ+fjZwa19tcnA3szamkqmU2iq4i4AZwN/k/DN465AMpBkxE4DHSQH9qohYlMfUvwDcDSwGboqIJ/IxXwO+ImkZaQz+ir4a4zF3M2tfdXxCtbdVcCPinB7yNpGmQ/Z0njuAO3rIf5Y0Vl+Kg7uZtS8vP2BmVjH1f4hpSHFwN7O25Zd1mJlVUcl1Y1qRg7uZtS25525mVjHlH1BqSaXmuUsaXSbPzKx1KN1QLZNaUNmHmB4omWdm1jrqt/zAkNPnsIykd5JWHhsr6QjefFRrIjCuwW0zM2usrsFuQOP0N+b+UeAc0iI1PyjkbwS+0aA2mZk1XjvPc4+IucBcSZ+KiJub1CYzs6Zo+9kyEXGzpI8DhwJjCvnfaVTDzMwarsLBvexsmcuA00mrnom04M2+DWyXmZnthLKzZd4fEWcDL0fEt4FjeOuC8mZmLUdRLrWisg8xvZa/vippb2AtsF+9GzNMYpdhw+t92rfZd8TGhtdRs8uwLU2pZ3gT/758qXNC0+oaOXZ7/zvVSVf5dbt32qZtY/rfqQ62dzX+/1PN6u27NK2uugi8/ABwu6RJwN8BD5O+LT9uWKvMzJqhRXvlZZS9ofqX+ePNkm4HxkTEhsY1y8ys8Vp1yKWM/h5iOq2PMiLilvo3ycysSdo1uAP/q9t27Vuh/NnB3cxaV7sG94g4F0DSGOBTwPTCMRX+tphZ1bXyTJgyyt5Q/VdgPelmam36R4W/LWbWFjxbhmkRcfKOVCBpOLAAeD4iPrEj5zAza4Qq99zLPsT0C0nv2cE6vggs3sFjzcwap8JL/pYN7h8EHpL0lKRFkh6TtKi/gyRNAz6O58Sb2VBT8unUVu3dlx2W+dgOnv/vga8CvT66JmkOMAdg2tTmPU1nZtaqvfIyyj7E9KuBnljSJ4CXIuIhScf3ce7LgcsBZhw+qsLfajMbalThl3WUHZbZER8APilpOXADcIKkaxtYn5mZZQ0L7hHx9YiYFhHTgTOAn0TE5xpVn5nZgFX4hmrZMXczs2pp4ZulZTQluEfEfcB9zajLzKw0B3czswpycDczqxbh2TJmZtVT54eYJE2SNE/SEkmLJR0j6UZJC3NaLmlh3vd3C/kLJXVJmpHLjsoPii6TdIkk5fwpku6RtDR/ndxXexzczax91Xe2zMXAXRFxEHA4sDgiTo+IGRExA7iZvEx6RFxXyD8LWB4RC/N5LiU92HlATrV1vS4E7o2IA4B783avHNzNrH3VKbhLmggcB1wBEBFbI2J9oVzAZ4Hrezj8zFq+pL2AiRHxQEQEcA1wat7vFGBu/jy3kN+jITXmPoJhTB42ruH1DH/jfd+Nt7qrWb8/mzd4+GqMbl5dXc2r64Utk5pW15bO5vzXa1Y9AHevPKRpdcGddTnLAKZC7iZpQWH78vx0fc3+wGrgKkmHAw8BX4yIzbn8WGBVRCzt4dynkwI3wFSgo1DWkfMA9oyIlQARsVLSHn012D13M2tf5XvuayJiZiFd3u1MI4AjgUsj4ghgM28dNnmjd14k6X3AqxHxeC2rl1YOmIO7mbWnSLNlyqQSOoCOiJift+eRgj2SRgCnATf2cNwZvDXodwDTCtvTgBfy51V52KY2fPNSXw1ycDez9lWnMfeIeBFYIenAnHUi8GT+fBKwJCKKwy1IGgZ8hrT2Vu08K4GNkmblcfqzgVtz8W3A7Px5diG/R0NqzN3MrJnqvPzABcB1kkYBzwLn5vzuvfOa40i9/We75Z8PXA2MJd1cqN1guAi4SdJ5wHOkXwy9cnA3s/ZVx+CepzLO7CH/nF72vw+Y1UP+AuCwHvLXkv4iKMXB3czaUwuv+FiGg7uZtSXhVSHNzCrJwd3MrIoc3M3MKsjB3cysYvwmJjOzinJwNzOrniq/rMPB3czalodlzMyqxg8xmZlVlIO7mVm1+AlVM7OKUld1o7uDu5m1J4+5m5lVk4dlzMyqyMG9Wl7obN6TC/dvPrgp9WzoHNuUegBe2DKpaXU9vWH3ptXV1eO7iRvjta0jm1LPS89Pbko9ACPXtl44cc/dzKyKHNzNzComvPyAmVnleJ67mVlVRXWju4O7mbUt99zNzKqm4g8xDWvUiSXtI+mnkhZLekLSFxtVl5nZjlBXudSKGtlz3w78SUQ8LGkX4CFJ90TEkw2s08ystFYN3GU0LLhHxEpgZf68UdJiYCrg4G5mgy/wDdWdJWk6cAQwv4eyOcAcgHdN9S0AM2ueKt9QbdiYe42kCcDNwJci4pXu5RFxeUTMjIiZu+86vNHNMTN7U5RMLaihXWVJI0mB/bqIuKWRdZmZDYQfYtpBkgRcASyOiB80qh4zsx0SUemXdTRyWOYDwFnACZIW5vRbDazPzGxgPCwzcBFxPzRxDVUzswGq8rBMw2+ompkNSQF0RblUgqRJkuZJWpIf3jxG0o2FkYvlkhYW9n+vpAfyQ56PSRqT84/K28skXZKHuJE0RdI9kpbmr30u1u/gbmbtq77DMhcDd0XEQcDhpPuNp0fEjIiYQZpccguApBHAtcDnI+JQ4HhgWz7PpaTp4QfkdHLOvxC4NyIOAO7N271ycDeztqUol/o9jzQROI40iYSI2BoR6wvlAj4LXJ+zPgIsiohH8/5rI6JT0l7AxIh4ICICuAY4NR9zCjA3f55byO+Rg7uZtS11RakE7CZpQSHN6Xaq/YHVwFWSHpH0Y0njC+XHAqsiYmnefjcQku6W9LCkr+b8qUBH4biOnAewZ37yv7YCwB59XZsfCTWz9jSwIZc1ETGzj/IRwJHABRExX9LFpGGTb+byM3mz117b/4PAbwCvAvdKegh424OeA2plgXvuZtaW0kNMUSqV0AF0RERtiZV5pGBfG18/Dbix2/4/i4g1EfEqcEfevwOYVthvGvBC/rwqD9uQv77UV4OGVM99bdcIrts4peH1PPbq4Q2vo2bdtvH971QHj6/dqyn1AKzZ0JxrAti+pXk/onplZNPqGr22Of2qvZ9q3rKH0cQlFp+p14nq1OSIeFHSCkkHRsRTwIm8uUjiScCSiCgOt9wNfFXSOGAr8CHghxGxUtJGSbNIa3GdDfzffMxtwGzgovz11r7aNKSCu5lZM5XslZd1AXCdpFHAs8C5Of8M3jokQ0S8LOkHwIOkYZc7IuLfc/H5wNXAWODOnCAF9ZsknQc8B3ymr8Y4uJtZe6rz06cRsRB427h8RJzTy/7XkqZDds9fABzWQ/5a0l8EpTi4m1mbqvbaMg7uZta+/LIOM7OKCb9mz8ysmtxzNzOroOrGdgd3M2tf6qruuIyDu5m1p6BuDzENRQ7uZtaWROmlBVqSg7uZtS8HdzOzCnJwNzOrGI+5m5lVk2fLmJlVTnhYxsyscgIHdzOzSqruqIyDu5m1L89zNzOrIgd3M7OKiYDO6o7LOLibWftyz7051m0dz7UvzGp4Pb9cvWvD66jZunZscyrarubUAwzb2ry6xq1pXl1TlnQ2ra4Jy9Y1p6ImzuPWlm1Nq6tuHNzNzComAL9D1cysagLCY+5mZtUS+IaqmVkleczdzKyCHNzNzKqm2guHDWvkySWdLOkpScskXdjIuszMBiRIU0XLpBbUsOAuaTjwI+BjwCHAmZIOaVR9ZmYDFlEutaBGDsscDSyLiGcBJN0AnAI82cA6zcxK8vIDO2oqsKKw3QG8r/tOkuYAcwDG7LlLA5tjZlYQEBWe597IMfeenht/2983EXF5RMyMiJmj3tGkR/XNzCA9oVomtaBG9tw7gH0K29OAFxpYn5nZwLToeHoZjey5PwgcIGk/SaOAM4DbGlifmVl5EXWdLSNpkqR5kpZIWizpGEk3SlqY03JJC/O+0yW9Vii7rHCeoyQ9lmcZXiJJOX+KpHskLc1fJ/fVnoYF94jYDnwBuBtYDNwUEU80qj4zswGr72yZi4G7IuIg4HBgcUScHhEzImIGcDNwS2H/Z2plEfH5Qv6lpPuQB+R0cs6/ELg3Ig4A7s3bvWroQ0wRcQdwRyPrMDPbMUF01meZZ0kTgeOAcwAiYiuwtVAu4LPACf2cZy9gYkQ8kLevAU4F7iTNNjw+7zoXuA/4Wm/nauhDTGZmQ1Ztyd9yN1R3k7SgkOZ0O9v+wGrgKkmPSPqxpPGF8mOBVRGxtJC3X973Z5KOzXlTSfcrazpyHsCeEbESIH/do6/L8/IDZta+yk+FXBMRM/soHwEcCVwQEfMlXUwaNvlmLuiz5zYAAASJSURBVD8TuL6w/0rgXRGxVtJRwL9KOpSSswzLcM/dzNpSANEVpVIJHUBHRMzP2/NIwR5JI4DTgBvfqDvi9YhYmz8/BDwDvDufZ1rhvMVZhqvysE1t+Oalvhrk4G5m7SnyyzrKpH5PFS8CKyQdmLNO5M2n8U8ClkTEG8MtknbPS7QgaX/SjdNn83DLRkmz8jj92cCt+bDbgNn58+xCfo88LGNmbateN1SzC4Dr8tTvZ4Fzc/4ZvHVIBtLN1+9I2g50Ap+PiNqLdc8HrgbGkm6k3pnzLwJuknQe8Bzwmb4aoxhCk/glrQZ+NcDDdgPWNKA5g62K11XFa4JqXtdQv6Z9I2L3nTmBpLtI11nGmog4uf/dho4hFdx3hKQF/dzoaElVvK4qXhNU87qqeE3txmPuZmYV5OBuZlZBVQjulw92AxqkitdVxWuCal5XFa+prbT8mLuZmb1dFXruZmbWjYO7mVkFtWxwl3SypKfymsd9Ln3ZKiTtI+mneS3oJyR9cbDbVC+ShudFkm4f7LbUS0/rdw92m+pB0pfzz9/jkq6XNGaw22QD15LBPT+2+yPgY8AhwJmSDhncVtXFduBPIuJgYBbwRxW5LoAvktb1r5K3rd89yO3ZaZKmAn8MzIyIw4DhpCcsrcW0ZHAHjgaWRcSzed3kG0hrHbe0iFgZEQ/nzxtJwWJq30cNfZKmAR8HfjzYbamXwvrdV0Bavzsi1g9uq+pmBDA2L3g1Dr8esyW1anCfCqwobBfXPK4ESdOBI4D5fe/ZEv4e+CpQpVfN97d+d0uKiOeB75HWLlkJbIiI/xjcVtmOaNXgXrc1j4ciSRNIr+T6UkS8Mtjt2RmSPgG8lJc1rZLa+t2XRsQRwGb6ee1ZK8jv5TwF2A/YGxgv6XOD2yrbEa0a3DuAfQrbxTWPW5qkkaTAfl1E3NLf/i3gA8AnJS0nDZ+dIOnawW1SXfS6fneLOwn4ZUSsjohtpHd+vn+Q22Q7oFWD+4PAAZL2y8trnkFa67il5fWbryC9WPcHg92eeoiIr0fEtIiYTvp3+klEtHxPsJ/1u1vZc8AsSePyz+OJVOBGcTtqyfXcI2K7pC8Ad5Pu5l8ZEU8McrPq4QPAWcBjkhbmvG/kF43b0NPb+t0tK78ibh7wMGn21iN4KYKW5OUHzMwqqFWHZczMrA8O7mZmFeTgbmZWQQ7uZmYV5OBuZlZBDu425Ei6WtKnB7sdZq3Mwd3MrIIc3G3QSJqe10H/x7x++H9IGtttnxPzwlyPSbpS0uicv1zStyU9nMsOGpyrMBuaHNxtsB0A/CgiDgXWA5+qFeSXRFwNnB4R7yE9UX1+4dg1EXEkcCnwp01rsVkLcHC3wfbLiKgttfAQML1QdmAufzpvzyWtoV5zSy/HmbU9B3cbbK8XPnfy1vWOelrauadjux9n1vYc3G0oWwJMl/Trefss4GeD2B6zluHgbkNWRGwhrbT4z5IeI73J6bLBbZVZa/CqkGZmFeSeu5lZBTm4m5lVkIO7mVkFObibmVWQg7uZWQU5uJuZVZCDu5lZBf1/N9n17Hr+PXoAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"h_prime.isel(time=0).plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Integrate O2 between `h_prime` and `h`"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"def interval_overlap(interval_1, interval_2):\n",
" \"\"\"Compute the overlap between two intervals.\"\"\" \n",
" overlap = (np.min([np.max(interval_1), np.max(interval_2)]) - \n",
" np.max([np.min(interval_1), np.min(interval_2)]))\n",
" return max([0.0, overlap])\n",
"\n",
"def compute_vertical_weights(coord_interval, coord_bounds):\n",
" \"\"\"Compute averaging weights for remapping a peicewise-constant field on a\n",
" 1D coordinate to a single interval.\n",
" \n",
" \n",
" Example: \n",
" [1]: weights = compute_vertical_weights(coord_interval=np.array([0, 20]), \n",
" coord_bounds=np.array([[0, 10], \n",
" [10, 20], \n",
" [20, 30]]))\n",
"\n",
" [2]: weights\n",
" \n",
" array([0.5, 0.5, 0. ])\n",
" \n",
" \n",
" Parameters\n",
" ---------\n",
" coord_interval : numpy.array\n",
" A two-element vector with the coordinate values over which to compute weights. \n",
" For example\n",
" \n",
" coord_bounds : numpy.array\n",
" A N x 2 element vector with the coordinate bounds.\n",
" \"\"\"\n",
" \n",
" weights = np.zeros(coord_bounds.shape[0])\n",
" for i in range(coord_bounds.shape[0]):\n",
" weights[i] = interval_overlap(coord_interval, coord_bounds[i, :])\n",
" return weights / weights.sum()\n"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 2min, sys: 10.4 s, total: 2min 10s\n",
"Wall time: 1min 59s\n"
]
}
],
"source": [
"%%time\n",
"nmolcm2_to_molm2 = 1e-9 * 1e4\n",
"\n",
"o2_heave = xr.DataArray(np.ones((nl, nj, ni)) * np.nan, \n",
" dims=('time', 'nlat', 'nlon'), \n",
" attrs={'units': 'mol m$^{-2}$'})\n",
"\n",
"sigma_remap_weights = xr.full_like(dss.sigma, fill_value=0.)\n",
"sigma_remap_weights.name = 'weights'\n",
"\n",
"o2_work_array = xr.full_like(dss.time, fill_value=0.)\n",
"o2_work_array.name = 'work'\n",
"\n",
"for j in range(len(dss.nlat)):\n",
" for i in range(len(dss.nlon)):\n",
" \n",
" if not np.isnan(sigma_z_clim[j, i]):\n",
" \n",
" # compute vertical distance between z=h and z=h_prime\n",
" dh = np.fabs(h - h_prime[:, j, i]) # cm\n",
"\n",
" # loop over time and compute mean O2 in sigma range\n",
" for l in range(len(dss.time)): \n",
" # compute weights \n",
" sigma_remap_weights[:] = compute_vertical_weights((sigma_z_t[l, j, i], sigma_z_clim[j, i]), sigma_bounds)\n",
" \n",
" # compute mean o2 between sigma levels\n",
" o2_work_array[l] = (sigma_remap_weights * dss.O2[l, :, j, i]).sum('sigma')\n",
"\n",
" # compute the vertical integral \n",
" o2_heave[:, j, i] = o2_work_array * dh * nmolcm2_to_molm2\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2ae080870eb8>"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAdqklEQVR4nO3dfbRdZWHn8e8v9+YFUEBBO5qgwQVV0da3DNhiu8ZSLU5t44y4jLaWupjJ6gutfbEt9kWFadcaOh2t07LsUMAipYNrQNtMG8VOqdN50ZTwUmNA2jRSuWCXhIS8kZDce3/zx9kXDseTc/a9d+997j3791lrL87ZZ5/neXYS7u8++9nPs2WbiIiIXitG3YCIiFiaEhAREdFXAiIiIvpKQERERF8JiIiI6Gty1A3oduZzJ7z+rJW11/OEZ2uvY85xN/NH/GRD9QAcmlndWF3HZiYaq2tmtrnfl2ZnGqprWs3UA2imsap48htTe2w/bzFl/MAbT/Fje8s1+q4vP3m77YsXU99ytKQCYv1ZK/nb219Uez33Hjtaex1zHp4+vZF6dh97fiP1AHzp8Zc0VtfXDz6nsbr2Hjq5sbqOHFjTSD0r9tb/C9ecVfubC6MHrvyFf1psGY/tnSn982biBf9w5mLrW46WVEBERDTFwCzNXU1YjhIQEdFKxhx3g9fFlqEERES0VnoQgyUgIqKVjJnJUkMDJSAiorVmSUAMkoCIiFYyMJOAGCgBERGtlR7EYAmIiGglA8czBjFQAiIiWsk4l5iGSEBERDsZZpIPAyUgIqKVOjOpY5AERES0lJihufWjlqMERES0UmeQOgExSAIiIlqpMw8iATFIAiIiWms2PYiBEhAR0UrpQQyXgIiIVjJiJk9dHigBERGtlUtMgyUgIqKVjDjm5p55vhwlICKilToT5XKJaZD86UREa80Uk+WGbWVIuljSA5J2Sbqiz+erJX2q+HybpPXF/vWSjki6t9j+oOs7r5O0o/jOf5HU6DWxJdWDmGGWA7NHaq/n0Zln117HnIeOn9FIPfcePKuRegAePPDcxuras/9ZjdV1fO+axuqaPNDM72aTh5r7ebLyUGNVVcIWM67m70HSBHAN8CZgCrhT0hbb93Uddhmwz/Y5kjYBVwPvLD77R9uv7lP0x4HNwJeArcDFwGcraXQJ6UFERGvNolJbCecDu2zvtn0MuAXY2HPMRuDG4vWtwEWDegSSXgCcavuLtg18EnjbfM9xMZZUDyIioimdQerSPwLPlLS96/21tq/ter8WeKjr/RRwQU8ZTx1je1rSfmDuEsPZku4BDgC/bvt/F8dP9ZS5tmyDq5CAiIhWmucg9R7bGwZ83q8n0LuY+ImO+QbwItuPSXod8KeSXlGyzFolICKitWaqmwcxBXQPBK4DHjnBMVOSJoHTgL3F5aMnAWzfJekfgW8vjl83pMxaZQwiIlppbiZ1ma2EO4FzJZ0taRWwCdjSc8wW4NLi9SXAHbYt6XnFIDeSXgKcC+y2/Q3goKTXF2MVPwb82eLPvLz0ICKitWYruoupGFO4HLgdmABusL1T0lXAdttbgOuBmyTtAvbSCRGA7wWukjQNzAA/YXtv8dlPAn8EnETn7qXG7mCCmgNC0s8D/47OdbMdwHttH62zzoiIMjqL9VV3EcX2Vjq3onbv+2DX66PAO/p87zbgthOUuR14ZWWNnKfaLjFJWgv8LLDB9ivppOqmwd+KiGiGEcc9UWprq7ovMU0CJ0k6DpxMwwMsEREnYlPZRLlxVdufju2Hgd8Bvk7nNq79tj/fe5ykzZK2S9q+57E8QjwimlJuklzJiXJjqc5LTM+hM3PwbOCFwCmSfrT3ONvX2t5ge8OZZyTNI6IZptODKLO1VZ1n/v3A12w/avs48Gngu2usLyJiXiq8zXUs1TkG8XXg9ZJOBo4AFwHbB38lIqIZRnlg0BC1BYTtbZJuBe4GpoF7gGsHfysiohkGjpdfi6mVav3Tsf0h4EN11hERsTDln/XQVonPiGglU91M6nGVgIiI1koPYrAERES0kq30IIZIQEREK3UGqdu7jEYZCYiIaKnqnkk9rpZUQDzhFfzdsZNqr2fH0bOGH1SRew68qJF6/n7v8xqpB+DRPac2VpceX9lYXav2N/fDYvJIQ/U80Uw90Nw5VaUzSJ0xiEGWVEBERDSpzbOky0hAREQrZSb1cAmIiGit2fQgBkpAREQr2XB8NgExSAIiIlqpc4kpATFIAiIiWiszqQdLQEREK+U21+ESEBHRUrnENEwCIiJaq83Pmy4jARERrdS5iylrMQ2SgIiIVspEueFyAS4iWmsWldrKkHSxpAck7ZJ0RZ/PV0v6VPH5Nknrez5/kaRDkt7fte9BSTsk3Stp+yJPd97Sg4iIVqryLiZJE8A1wJuAKeBOSVts39d12GXAPtvnSNoEXA28s+vzjwKf7VP8G23vqaSh85QeRES01qxXlNpKOB/YZXu37WPALcDGnmM2AjcWr28FLpIkAElvA3YDOys5sYokICKilWwx7RWlthLWAg91vZ8q9vU9xvY0sB84Q9IpwK8AV/ZrJvB5SXdJ2jzPU1y0XGKKiNaaxyWmM3vGAK61fW3X+34Fuef9iY65Evio7UNFh6LbhbYfkfR84C8lfdX235Rt9GIlICKileY5BrHH9oYBn08B3U8iWwc8coJjpiRNAqcBe4ELgEsk/TZwOjAr6ajt37f9CIDtb0r6DJ1LWQmIiIi6VXib653AuZLOBh4GNgHv7jlmC3Ap8EXgEuAO2wa+Z+4ASR8GDtn+/eLS0wrbB4vXbwauqqrBZSQgIqKVqpwHYXta0uXA7cAEcIPtnZKuArbb3gJcD9wkaRednsOmIcV+G/CZ4rLTJPAntj9XSYNLSkBERGtVudSG7a3A1p59H+x6fRR4x5AyPtz1ejfwqsoauAAJiIhoJRum88CggZZUQDw5u5J/OPZttdez+8jzaq9jzj/sa6aux/Y9q5F6ADiwsrGqJo409z/wiunGqmLySDP1rDzcTD0Ak0d6b9pZ+rLUxmBLKiAiIpqStZiGS0BERGs5ATFQAiIiWivPgxgsARERrWRnDGKYBEREtJSYyV1MAyUgIqK1MgYxWK3xKel0SbdK+qqk+yV9V531RUSUNbcWU5mtreruQXwM+JztSyStAk6uub6IiHLcGYeIE6stICSdCnwv8OMAxUM0jtVVX0TEfOUupsHqvMT0EuBR4BOS7pF0XbEi4TNI2ixpu6TtB/cdr7E5ERFPczFIXWZrqzrPfBJ4LfBx268BDgPf8iBv29fa3mB7w7Of09wSDhERdrmtreoMiClgyva24v2tdAIjImJJsFVqa6vaAsL2PwMPSXppsesi4L666ouImI9O7yABMUjddzH9DHBzcQfTbuC9NdcXEVFam29hLaPWgLB9LzDoOa4RESPT5vGFMjKTOiJayYjZFt+hVEYCIiJaKx2IwRIQEdFOHq+1mCQ9t8Rhs7YfL1tmAiIi2mu8uhCPFNug1JsAXlS2wARERLTWOPUggPuLScknJOme+RSYgIiIVjIwOztWAVFmtex5rai9pALi8Oxq/vbAS2qv5+/3P6/2Oubsf2JNI/XMPtncX+WK4839T6XZxqpi8nBzda082Ew9a/bNNFMRMPlEg39ZVTAwRj0I20erOKZb7vGKiNYal7WYJL1J0h9KenXxfnMV5S6pHkRERKOWwQ//kn6KzkoVv17czfTqKgpNDyIiWqrcOkxlB7IlXSzpAUm7JH3LytWSVkv6VPH5Nknrez5/kaRDkt5ftswuj9p+3Pb7gTcD/7JUo4dIQEREe7nkNoSkCeAa4C3AecC7JJ3Xc9hlwD7b5wAfBa7u+fyjwGfnWeacv3jqlOwrgE8Ob/VwCYiIaCeDZ1VqK+F8YJft3cXTM28BNvYcsxG4sXh9K3CRJAFIehudBU13zrPMzqnYf1aUc2bx/vfKNHqYBEREtJhKbpw59+TLYusdBF4LPNT1fqrY1/cY29PAfuCM4kmbvwJcuYAye90w5PN5ySB1RLRX+UHqPbYHrUzdr5vRW/qJjrkS+KjtQ0WHYj5llmnHgiUgIqK9qruLaQo4q+v9OjrLXvQ7ZkrSJHAasBe4ALhE0m8DpwOzko4Cd5Uos1el92UlICKinaqdKHcncK6ks4GHgU3Au3uO2QJcCnwRuAS4w7aB75k7QNKHgUO2f78IkWFl9koPIiKiClVNgrM9Lely4HY6C+LdYHunpKuA7ba3ANcDN0naRafnsGkhZQ5pygcWey7dEhAR0V4VrsVkeyuwtWffB7teHwXeMaSMDw8rc8j3v1L22DISEBHRWhqfmdRPkbQB+DXgxXR+xguw7e+cb1kJiIhop5KT4Jahm4FfAnYAi1pBsVRASFpt+8lh+yIilg+N1WquXR4txjwWrWwP4ovAa0vsi4hYPsazB/EhSdcBfwU89Uu87U/Pt6CBASHpX9CZuXeSpNfw9C1UpwInz7eyiIglZZk9wqKk9wIvA1by9BkaqDYggB8AfpzOBI2PdO0/CPzqfCuLiFgyxuyBQV1eZfs7qihoYEDYvhG4UdLbbd9WRYUREUvFON7FBHxJ0nm271tsQaXGIGzfJukHgVcAa7r2X7XYBkREjMx4BsQbgEslfY3OGES9t7lK+gM6Yw5vBK6jM038b+dbWURE1O7iqgoqexfTd9v+Tklftn2lpP/MAgY8IiKWknG8xGT7n6oqq2xAHCn++4SkFwKPAWdX1YinKpleyZcfe2HVxX6LvYeauwHryf2rG6lnxeGJRuoBmDjS3MDeqv2NVcWqA83VddLemUbqWXmwmXoAJp+YbqyuSphKl9oYR2UD4s8lnQ78J+BuOn+019XWqoiIJoxhD6JKZQep/0Px8jZJfw6ssd3g73YREdUbx0tMVRo2Ue7fDvhsQTPzIiKWjDEKCEkHeeYZqXg/dxfTqfMtc1gP4od63s9VPldxAiIilq8xCgjbz666zGET5d4LIGkN8HZgfdd3xuiPNiLaRh7fS0ySXsXTT6r7G9tfXkg5K0oe96d0ehPHgUNdW0TE8jWrctsyIul9dJb8fn6x3SzpZxZSVtm7mNbZXtDkC0kTwHbgYdtvXUgZERF1GNMexGXABbYPA0i6ms7q278334LK9iD+n6SFLv70PuD+BX43IqI+LrktLwK6J8DM8PRK3PNStgfxBuDH57u2h6R1wA8CvwX8wkIaGBFRi/Edg/gEsE3SZ+j8rH4bcMNCCiobEG9ZSOHA7wK/DJxwdF3SZmAzwKrnz/surIiIhRvDgLD9EUlfAC6kExCX2r53IWWVnSg377U9JL0V+KbtuyT9qwFlXwtcC/Csb3/BGP51RcRSpTF8YJCkDcCv8fRdp/9eUn2ruS7QhcAPS/rXdJYIP1XSH9v+0RrrjIhou5uBXwJ2sMhn5tUWELY/AHwAoOhBvD/hEBFLynhes3jU9pYqCqqzBxERsXSN7yD1hyRdB/wVnZuKABa0NFIjAWH7C8AXmqgrIqK0CgNC0sXAx4AJ4Drb/7Hn89XAJ4HX0XlkwjttPyjpfIpxWDqDyh+2/ZniOw8CB+ncqjpte0OJprwXeBmwkqcvMS1oaaT0ICKivSoKiGJC8DXAm4Ap4E5JW3qeC30ZsM/2OZI2AVcD7wS+AmywPS3pBcDfSfoftucesPFG23vm0ZxX2V7ovLVnKDtRLiJirIjOXUxlthLOB3bZ3m37GHALsLHnmI3AjcXrW4GL1Lm96ImuMFjD4mPrS5LOW2QZQAIiItrKTy/YN2wDzpS0vWvb3FPaWuChrvdTxb6+xxSBsB84A0DSBZJ20rnz6Ce6AsPA5yXd1afOE3kDcK+kByR9WdIOSQtarC+XmCKivcr/rr5nyPX/fktZ9JZ+wmNsbwNeIenlwI2SPmv7KHCh7UckPR/4S0lftf03Q9q6oHXz+klARER7VTdIPQWc1fV+HfDICY6ZkjQJnAbsfUZz7PslHQZeCWy3/Uix/5vF0hnnAwMDYiETm09kSQXE9MwK9uw/pfZ6jh9aVXsdc1YcnmiknpUHmrtaOHm4sapY3eCDbdfsa25a7arHp4cfVIHJg8caqQdgxdHm6qpKhbe53gmcK+ls4GFgE/DunmO2AJfSWVn1EuAO2y6+81AxSP1i4KXAg5JOAVbYPli8fjNw1QnPRbrb9msHNbLMMd2WVEBERDSqooAofrhfDtxO5zbXG2zvlHQVnZ7AFuB64CZJu+j0HDYVX38DcIWk43RuS/0p23skvQT4jCTo/Kz+E9ufG9CMlw8ZaxCdXktpCYiIaCdXuxaT7a3A1p59H+x6fRR4R5/v3QTc1Gf/buBV82jCy0ocMzP8kKclICKivcZoJnWVYw9zEhAR0VpjutRGZRIQEdFeCYiBEhAR0U7L83GijUpAREQriVxiGiYBERGtlYAYLAEREe2VgBgoARER7ZWAGCgBERHtNL5PlKtMAiIi2isBMVACIiJaq8qlNsZRAiIiWiuXmAZLQEREO2Wi3FAJiIhorwTEQAmIiGilzKQeLgEREa2l2STEIAmIiGinjEEMlYCIiNbKJabBEhAR0V4JiIGWVEB4Rhw/sLr2eiYOTdRex5yVB9VIPav2NVINAKv3N/d/1ZrHm5vJtHrf8cbqmjzwZCP1rHjiWCP1AHDkaHN1VSQ9iMGWVEBERDQqATFQAiIi2slZamOYBEREtFLmQQyXgIiI9nISYpAERES0VnoQg60YdQMiIkbC89hKkHSxpAck7ZJ0RZ/PV0v6VPH5Nknri/3nS7q32P5O0r8pW2bdagsISWdJ+mtJ90vaKel9ddUVEbEQmi23DS1HmgCuAd4CnAe8S9J5PYddBuyzfQ7wUeDqYv9XgA22Xw1cDPxXSZMly6xVnT2IaeAXbb8ceD3w002fXETEIFUFBHA+sMv2btvHgFuAjT3HbARuLF7fClwkSbafsD1d7F/D032WMmXWqraAsP0N23cXrw8C9wNr66ovImJeTGeQuswGZ0ra3rVt7iltLfBQ1/spvvXn3VPHFIGwHzgDQNIFknYCO4CfKD4vU2atGhmkLq61vQbY1uezzcBmgInnnt5EcyIigHkNUu+xvWFQUX329ZZ+wmNsbwNeIenlwI2SPluyzFrVPkgt6VnAbcDP2T7Q+7nta21vsL1h4tmn1N2ciIinVTdIPQWc1fV+HfDIiY6RNAmcBux9RnPs+4HDwCtLllmrWgNC0ko64XCz7U/XWVdExHzMTZQrs5VwJ3CupLMlrQI2AVt6jtkCXFq8vgS4w7aL70wCSHox8FLgwZJl1qq2S0ySBFwP3G/7I3XVExGxIHZlDwyyPS3pcuB2YAK4wfZOSVcB221vofPz8CZJu+j0HDYVX38DcIWk48As8FO29wD0K7OSBpdU5xjEhcB7gB2S7i32/artrTXWGRFRXoVX9IufbVt79n2w6/VR4B19vncTcFPZMptUW0DY/j/0H2SJiFgSMpN6sCy1ERHtZCDPpB4oARER7ZV8GCgBERGtlUtMgyUgIqK1qrqLaVwlICKineaxUmtbJSAiopU6E+WSEIMsrYCYEZOPT9RezcrDzd19u/JgM/WsPtDcP/Q1+5p7kO/qx483VtfKRw81VpcOPtFMRU8ea6YewE8caayuyuSZ1AMtrYCIiGhQehCDJSAiop0yBjFUAiIiWqq6tZjGVQIiItorl5gGSkBERDu59ONEWysBERHtlR7EQAmIiGiv5MNACYiIaC3N5hrTIAmIiGgnk4lyQyQgIqKVhDNRbogERES0VwJioARERLRXAmKgBEREtFPGIIZKQEREa+UupsESEBHRUs4lpiESEBHRTiYBMUQCIiLaK1eYBlox6gZERIyK7FJbqbKkiyU9IGmXpCv6fL5a0qeKz7dJWl/sf5OkuyTtKP77fV3f+UJR5r3F9vyKTr2U9CAior0qusQkaQK4BngTMAXcKWmL7fu6DrsM2Gf7HEmbgKuBdwJ7gB+y/YikVwK3A2u7vvcjtrdX0tB5Sg8iItrJhpnZcttw5wO7bO+2fQy4BdjYc8xG4Mbi9a3ARZJk+x7bjxT7dwJrJK2u4AwXLQEREe1ll9vgTEnbu7bNPSWtBR7qej/FM3sBzzjG9jSwHzij55i3A/fYfrJr3yeKy0u/IUmLPON5WVKXmDQDKw/Vf/6rDtReRVddzdwlsWZfc6Ntq/ZPN1bX5L6jjdWl/Ycbq8sHDzVTz7FjjdQDMHvkSGN1Vab8JaY9tjcM+LzfD67ewgceI+kVdC47vbnr8x+x/bCkZwO3Ae8BPlmuyYuXHkREtJOBWZfbhpsCzup6vw545ETHSJoETgP2Fu/XAZ8Bfsz2Pz7VRPvh4r8HgT+hcymrMQmIiGgpg2fLbcPdCZwr6WxJq4BNwJaeY7YAlxavLwHusG1JpwN/AXzA9v+dO1jSpKQzi9crgbcCX1nUKc/TkrrEFBHRGFN2AHp4Ufa0pMvp3IE0Adxge6ekq4DttrcA1wM3SdpFp+ewqfj65cA5wG9I+o1i35uBw8DtRThMAP8T+MNKGlxSAiIi2qvCmdS2twJbe/Z9sOv1UeAdfb73m8BvnqDY11XWwAVIQEREe2WpjYESEBHRUlmsb5haB6mHTT2PiBgZA7Oz5baWqi0guqaevwU4D3iXpPPqqi8iYt7KT5RrpTovMT019RxA0tzU8/sGfisiohGu7C6mcVVnQPSben5B70HFlPXNAJOnPafG5kREdDG43ByH1qpzDKLM1HNsX2t7g+0NEyefUmNzIiJ6VDeTeizV2YMoM/U8ImJ0Wjy+UEadAfHU1HPgYTqzBt9dY30REeXZrb5DqYzaAuJEU8/rqi8iYt7Sgxio1oly/aaeR0QsDcYzM6NuxJKWmdQR0U5zy33HCSUgIqK9cpvrQAmIiGglA04PYqAERES0k50exBAJiIhorQxSDyYvodu8JD0K/NM8v3YmsKeG5ozaOJ7XOJ4TjOd5LfVzerHt5y2mAEmfo3OeZeyxffFi6luOllRALISk7bY3jLodVRvH8xrHc4LxPK9xPKeYv1qfBxEREctXAiIiIvoah4C4dtQNqMk4ntc4nhOM53mN4znFPC37MYiIiKjHOPQgIiKiBgmIiIjoa9kGhKSLJT0gaZekK0bdnipIOkvSX0u6X9JOSe8bdZuqImlC0j2S/nzUbamKpNMl3Srpq8Xf2XeNuk1VkPTzxb+/r0j6b5LWjLpNMRrLMiAkTQDXAG8BzgPeJem80baqEtPAL9p+OfB64KfH5LwA3gfcP+pGVOxjwOdsvwx4FWNwfpLWAj8LbLD9SjrPctk02lbFqCzLgADOB3bZ3m37GHALsHHEbVo029+wfXfx+iCdHzhrR9uqxZO0DvhB4LpRt6Uqkk4Fvhe4HsD2MduPj7ZVlZkETpI0CZxMHhXcWss1INYCD3W9n2IMfpB2k7QeeA2wbbQtqcTvAr8MjNPKaC8BHgU+UVw6u07SKaNu1GLZfhj4HeDrwDeA/bY/P9pWxags14BQn31jc7+upGcBtwE/Z/vAqNuzGJLeCnzT9l2jbkvFJoHXAh+3/RrgMLDsx8IkPYdOb/xs4IXAKZJ+dLStilFZrgExBZzV9X4dY9INlrSSTjjcbPvTo25PBS4EfljSg3QuBX6fpD8ebZMqMQVM2Z7r4d1KJzCWu+8Hvmb7UdvHgU8D3z3iNsWILNeAuBM4V9LZklbRGUTbMuI2LZok0bmmfb/tj4y6PVWw/QHb62yvp/P3dIftZf8bqe1/Bh6S9NJi10XAfSNsUlW+Drxe0snFv8eLGIPB91iYZfk8CNvTki4Hbqdzl8UNtneOuFlVuBB4D7BD0r3Fvl+1vXWEbYoT+xng5uKXlN3Ae0fcnkWzvU3SrcDddO6qu4csu9FaWWojIiL6Wq6XmCIiomYJiIiI6CsBERERfSUgIiKirwRERET0lYCIJUfSH0m6ZNTtiGi7BERERPSVgIiRkbS+eI7CHxbPH/i8pJN6jrmoWAxvh6QbJK0u9j8o6UpJdxefvWw0ZxExvhIQMWrnAtfYfgXwOPD2uQ+KB9X8EfBO299BZ+b/T3Z9d4/t1wIfB97fWIsjWiIBEaP2Ndtzy4rcBazv+uylxed/X7y/kc4zGOZ8+gTfi4gKJCBi1J7sej3DM9cH67ese7/v9n4vIiqQgIil7KvAeknnFO/fA/yvEbYnolUSELFk2T5KZ4XU/y5pB50n0v3BaFsV0R5ZzTUiIvpKDyIiIvpKQERERF8JiIiI6CsBERERfSUgIiKirwRERET0lYCIiIi+/j+y87Seb2IW9gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"o2_heave.std('time').plot()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2ae0674a4668>"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAco0lEQVR4nO3deZQlZZ3m8e9TC6siCNq2gICOgqIiWK6ojaIt7u1yjqJoS9tN24MK7uL0SDtMT7u3jK3SBQ24MOWMgIq44DICPW2DFkuzj+2hFEtwigKhiiqoJfOZPyJS0zTz3sjMiLiRGc/nnDjkvTduvL97i8qn3njfeEO2iYiImGrJqAuIiIhuSkBERMS0EhARETGtBEREREwrAREREdNaNuoCJtvzgUu8776dKmneNOoCGqAWP1Wb31+bnyvm54prtqy3/aD5HOP5z97Vd9w5VrW9i2wfNZ/2FqJO/Tbed99lfPdb8/ozr2RJi78IlrbUVqufSe21tYylrbW1XJ366xADLP3Df//5fI9xx51j/Oiih1Vtb6/5trcQ5W9ERPSSgXHGR11GpyUgIqKXjNnmaqeY+ioBERG9lR7EYAmIiOglY8ay1NBACYiI6K1xEhCD5DqIiOglA2O40jaMpDMlrZN03QyvP0DS1yX9m6TrJR1b9+dpQgIiInprHFfaKjgbGHSdxPHADbYPAY4APi5ph3l/gIblFFNE9JKBbTWNQdi+VNL+Q5q7vyQB9wPuBLbX0niDEhAR0UuuePqotJek1ZMer7S9chbN/QNwAXArcH/g1bY7P4UqARER/WQYq96BWG97xTxaez5wNfAc4BHAdyX9s+0N8zhm4zIGERG9VFxJXW2rwbHA+S78FFgDHFTPoZuTgIiInhJjFbca3AIcCSDpD4ADgZvrOHCTcoopInqpGKSuZ+FJSasoZiftJWktcDKwHMD2acApwNmSrqVYpPi9ttfX0niDEhAR0UvFdRD1BITto4e8fivwx7U01qIERET01nhNPYjFKgEREb1UZw9isUpAREQvGTGWeToDJSAiordyimmwBERE9JIRW93eLW0XogRERPRScaFcTjENkoCIiN7KIPVgnQoIQyt3eJrFAl3ztqW1ltqTu3DN3yzWAJqXbS3+AqzrorO22GLM6UEM0qmAiIho03h6EAMlICKil4pB6vwKHCTfTkT0Ugaph0tARERvjS2wcZO2JSAiopdyJfVwCYiI6K3xzGIaqNFvR9LbJV0v6TpJqyTt1GR7ERFVFYv1Lam09VVjn1zS3sDbgBW2HwssBV7TVHsREbNhxDYvrbT1VdOnmJYBO0vaBuwC3NpwexERldjkQrkhGvt2bP8S+BjFvVhvA+62/Z2p+0k6TtJqSavvuKOm24NHRAwlxitufdXkKaY9gJcBBwAPBXaVdMzU/WyvtL3C9oo990yaR0Q7iqV9llTa+qrJT/5cYI3t221vA84Hnt5gexERs5JB6sGaHIO4BXiqpF2Ae4EjgdUNthcRUZlRbhg0RGMBYftySecCVwLbgauAlU21FxExGwa2ZS2mgRr9dmyfDJzcZBsREXOj3A9iiMRnRPSSyZXUwyQgIqK30oMYLAEREb1kKz2IIRIQEdFLxSB1f5fRqCIBERE9lXtSD9PLgNhGS3eMB7a1dXP6Fudz39fi1MCttPcvvG0t/rK4z8tbaWfT+I6ttAPtfabC/Jd1KwapMwYxSC8DIiIC6PVV0lUkICKil3Il9XAJiIjorfH0IAZKQEREL9mwbTwBMUgCIiJ6qTjFlIAYJAEREb2VK6kHS0BERC9lmutwCYiI6KmcYhomARERvdXn+01XkYCIiF4qZjF1ay0mSZ+CmZd6sP22FstJQEREP9V5oZykM4EXA+tsP3aGfY4APgksB9bb/qNpduvUbZkTEBHRWzWeYjob+Afg89O9KGl34DPAUbZvkfTg6faz/bkp77t/8bTvqavQ2cgITUT00sQspirb0GPZlwJ3DtjltcD5tm8p91836HiSHivpKuA64AZJV0g6uPKHq0kCIiJ6a9xLKm3AXpJWT9qOm2VTjwL2kHRx+cv+DUP2Xwm8w/Z+th8GvBM4ffafcH5yiikieskW26tPc11ve8U8mlsGPBE4EtgZ+FdJl9n+yQz772r7B7+t1RdL2nUe7c9JAiIieqvFC+XWUoTMJmCTpEuBQ4CZAuJmSf8Z+EL5+BhgTfNl/q6cYoqIXqpzDKKCrwHPlLRM0i7AU4AbB+z/Z8CDgPOBr5Q/H1tHIbORHkRE9FaN01xXAUdQjFWsBU6mmM6K7dNs3yjp28A1wDhwhu3rZjqe7V8DrV7zMJ0ERET0Up3XQdg+usI+HwU+WuV4klYA7wf2Z9LvaduPn2OJc5KAiIje6vBSG+cA7waupehxjEQCIiJ6yYbt3b1h0O22Lxh1EZ0KCGO2zbwMSW02j7f3r4ZNbucr3jS+QyvtAGz0Tq21tWl8x9ba2ji2c3ttjbfzHd4z1t6f1d3b2/v+4MpajtLh5b5PlnQG8H1gy8STts9vs4hOBURERFvqHINowLHAQRQD3ROnmEwxq6k1CYiI6C13NyAOsf24URfR2RNwERFNG0eVthG4TNJjRtHwZOlBREQv2Z0eg3gG8KeS1lCMQYhiVddMc42IaJ4Y6+4spqNGXQAkICKix7o6BmH756OuARoeg5C0u6RzJd0k6UZJT2uyvYiIqlpei6kSSUPn71bZpy5N9yBOBb5t+1WSdgB2abi9iIhqXIxDdMyjJV0z4HUBD2irmMYCQtJuwLOANwLY3gpsbaq9iIjZ6uBSGwdV2Ges8SpKTfYgHg7cDpwl6RDgCuCEcj303yjvzHQcwN57d3bAKCIWGXdwkLorYw8Tmvx2lgGHAZ+1fSiwCXjf1J1sr7S9wvaKB+7ZrT+siFjc7GpbXzX5G3ktsNb25eXjcykCIyKiE2xV2vqqsYCw/SvgF5IOLJ86ErihqfYiImaj6B0kIAZpehbTW4FzyhlMNzOCW+ZFRMykq1dSS3oF8GHgwRQzlyaupN6tzToaDQjbVwMrmmwjImKuOjy+8BHgJbYH3be6cbmSOiJ6yYjxjs1imuT/jTocIAERET3W3Q4EqyX9T+Cr5IZBEREtc3fXYgJ2AzYDfzzpudwwKCKiNR3tQtjuxISezp6Ai4hoWlenuUp6lKTvS7qufPx4SX/ddh0JiIjoJQPj46q0jcDpwEnANgDb1wCvabuITp1i2uYl/Gpsx8bb2TjefBsT7hprZwHbu8Z2baUdgDtbbOvu7e0tALxh+06ttXXXtpb+v9i6cyvtAGzY2t7fq1oY6O4YxC62fyT9Tn3b2y6iUwEREdGmDl8HsV7SIyhHSSS9Crit7SISEBHRX90NiOOBlcBBkn4JrAGOabuIBERE9FR311myfTPwXEm7AktsbxxFHQmIiOivjvYgJL1jymOAu4EryiWMWpGAiIh+Mng0M5SqWFFuXy8fvwj4MfBmSV+2/ZE2ikhARESPdTYg9gQOs30PgKSTKe6p8yyKu3MmICIiGtXRU0zAw4Ctkx5vA/azfa+kLTO8p3YJiIjor+4GxP8ALpP0tfLxS4BV5aB1azdeS0BERD91+EI526dI+ibwDIrzYG+2vbp8+XVt1ZGAiIje6uqFcpI+Bpxl+9RR1pGAiIj+6u4sppuAlZKWAWcBq2zfPegNkh5Y4bjjtu+qWkQCIiJ6Sx3tQdg+AzhD0oHAscA1kv4FON32D2Z4263lNij1llIMgFeSgIiIfjJdHqRG0lLgoHJbD/wb8A5Jf2l7upVdb7R96JBjXjWbGiot9y3p95ZpnO65iIiFQ8UgdZVt2JGkMyWtm7h/w4D9niRprFx8b9B+n6A4zfRC4L/ZfqLtD9t+CTBTCDxtaKHV9vmNqveD+NeKz0VELByuuA13NnDUoB3KHsGHgYsqHO864BDbf2n7R1Nee/J0b7B937CDVtlnsoEBIekhkp4I7CzpUEmHldsRQHsL9UdENGG84jaE7UuBO4fs9lbgPGBdhcpke/NvHkhLy6upmW6wWtLzJJ0u6Qnl4+MqtDHUsDGI5wNvBPYBPjHp+Y3A++soICJiJFq8DkLS3sDLgecAT6rwliMlvRJ4E8WyG2cBlwzY/z9SDGb/dTmb6Qnzq7gwMCBsfw74nKRX2j6vjgYjIrpiFrOY9pK0etLjlbZXzqKpTwLvtT025S5x07L9WkmvBq4FNgNH2/6XAW+5vZy++i5JH6JaCA1VaRaT7fMkvQg4GNhp0vP/pY4iIiJGonpArLe9Yh4trQC+VIbDXsALJW23/dXpdpb0SOAEilNSjwZeL+mqyaedpvjGxA+23yfprfOo9TcqBYSk0yjGHJ4NnAG8Cpg6cBIREdOwfcDEz5LOBi6cKRxKXwfeYvt7KlLlHRTLfR88w/G/Vh57L9vrbX+qjrqrzmJ6uu03AL+2/UGKqVL71lFARMSoyNW2oceRVlHM7DxQ0lpJb5L0ZklvnmNpT7b9PQAXPg78SYX3nTnH9qZV9UK5e8v/bpb0UOAO4IAB+8/JfV7OTVseUvdhf8+dY/drvI0J67fdv5V27ti2ayvtAPx6a3sT2O7aunNrbW3YstPwnepq69522rr3vuWttAOwvcW2amFqW2rD9tGz2PeNg16X9GDgeEkHU1R5A/AZ2/9e4fC1jrpX7UFcKGl34KPAlcDPgC/VWUhEROvquw6iFpIOpziVBPB54Ivlz5eXrw1Ta7VVB6lPKX88T9KFwE7DFo6KiOi6Dq7F9HHgT2xPXhLja5K+Avwj8JQh76+1BzEwICS9YsBr2D6/zmIiIlrVvYDYbUo4AGD7aklVzlefVGcxw3oQL5nyeOLrVPlzAiIiFq7uBYQk7WH711OefCAVhgRsXydpBfCfgP0ofsereMmPn20xwy6UO7YsbifglcD+k97Tva82IqKiqjOUWvb3wHckvYtivBfgiRRrOP19xWOcA7yb4iK7CguFzKzqLKavAndRFDyx2FP3vtqIiNno2A2DbK+UdCtwCsU1DxOzmP6r7a9XPMztti+oo56qAbGP7YErFc6kXMFwNfBL2y+eyzEiIprQwR4Eti8ELpzHIU6WdAbwfWDLpOPOekigakD8UNLjbF872wYoLhe/EdhtDu+NiGhOBwOiBsdS3GRoOb89xTSnMeOqAfEM4I2S1lAkUqVBD0n7AC8C/pbiUvGIiG7o5hhEHQ6x/bg6DlQ1IF4wx+N/EngPMOP0rHLd8uMA9nzoDnNsJiJiDhZnQFwm6TG2b5jvgapeKPfz2R5Y0ouBdbavKG8wNNOxVwIrAQ543P0W5x9XRHSS5jXHp36SBp5psf2JQa+XngH86WzP+Eynag9iLg4HXirphRRLhO8m6Yu2j2mwzYiIhWzibMuBFPd0mJiN9BLg0orHmNOEouk0FhC2T6K8qq/sQbwr4RARndKxcxblatlI+g5wmO2N5eO/Ab5c8RizPuMzkyZ7EBER3dXtQeqHAVsnPd5KcaFyq1oJCNsXAxe30VZERGXdDYgvAD8qF+kzxf2sP992EelBRER/dTQgbP+tpG8BzyyfOna6RfyaloCIiF4S3ZvFNMUuwAbbZ0l6kKQDbK9ps4AERET0U4fHICSdDKygmM10FsVV0V+kmB3amqp3lIuIWHw6dke5SV4OvBTYBGD7VgZccNyU9CAior862oMAttq2VPRxJLV30/lJOhUQm8Z25EcbH9F4O+u3tvddr7+vnbbuunfnVtoBuOe+HVtra8u9y1tra+ze9v46LNm8tJV2lm1q7yTBTve21lRtunqKCfhfkv4R2F3SXwB/BpzedhGdCoiIiFZ1NCBsf0zS84ANFOMQH7D93bbrSEBERD+5u7OYJL0d+PIoQmGyDFJHRH91d5B6N+AiSf8s6XhJfzCKIhIQEdFbE/elHra1zfYHbR8MHA88FLhE0vfariOnmCKivzo6BjHJOuBXwB3Ag9tuPD2IiOinqqeXRhAikv5K0sUU95XeC/iLudzPYb7Sg4iIXhKdnua6H3Ci7atHWUQCIiJ6q2sBIemB5Y8fmfIYANt3tllPAiIi+qtjAQFcwW+r0pTXDDy8zWISEBHRXx0LCNsHjLqGyRIQEdFPHV7NFUDSS4FnlQ8vtn1h2zVkFlNE9Fd3ZzF9CDgBuKHcTpD0d23XkR5ERPRWV5faAF4IPMH2OICkzwFXASe1WUR6EBHRW129krq0+6SfHzCKAtKDiIh+Gt06S1X8HXCVpB9QzGZ6Fi33HiABERF91tGAsL2qvJL6SRQB8V7bv2q7jpxiiohemriSuounmCQdDmywfQHFrUbfI2m/tutIQEREb2nclbYR+CywWdIhwLuBnwOfb7uIBERE9FONi/VJOlPSOknXzfD66yRdU24/LH/xD7LdtoGXAf/d9qkUPYlWJSAiordqPMV0NnDUgNfXAH9Ursh6CrByyPE2SjoJOAb4hqSlQHs3aC8lICKiv2rqQdi+FJhxIT3bP7T96/LhZcA+Qw75amAL8KZycHpv4KPDK6lXp2Yxbdq+A5eva34cZsPmnRpvY8KWze2Evje390e55N6lrbW17J6p65U1Z6dNrTXF8pbaWn5Pe+fPl2/u6JSgAWYxAL2XpNWTHq+0PawXMJM3Ad8ass9G4FTbY5IeBRwErJpje3PWqYCIiGhV9YBYb3vFfJuT9GyKgHjGkF0vBZ4paQ+KmwatpuhVvG6+NcxGTjFFRD+5WGqjylYHSY8HzgBeZvuOYbvb3gy8AviU7ZcDB9dTSXUJiIjopTavg5D0MOB84PW2f1LtLXoaRY/hG+Vz7Z3bLeUUU0T0l+sZN5G0CjiCYqxiLXAy5awj26cBHwD2BD4jCYpprINOWZ1IsbTGV2xfL+nhwA9qKXYWEhAR0Vt1XSVt++ghr/858OezON4lwCWTHt8MvG3OBc5RAiIi+qmDi/VJ+qTtEyV9nWmqs/3SNutpLCAk7UtxafhDgHGKaWGnNtVeRMRsdfB+EF8o//uxkVZRarIHsR14p+0rJd0fuELSd23f0GCbERGVdS0gbF9R/vcSSQ8qf759VPU0NovJ9m22ryx/3gjcSHE1YETE6JlikLrK1hIV/kbSeuAm4CeSbpf0gdaKmKSVaa6S9gcOBS6f5rXjJK2WtHr73ZvbKCciAujkct8nAocDT7K9p+09gKcAh0t6e6uV0EJASLofcB5wou0NU1+3vdL2Ctsrlj1gl6bLiYj4rZrWYqrRG4Cjba/5TYnFDKZjytda1egsJknLKcLhHNvnN9lWRMRsTFwo1zHLba+f+qTt28vfp61qchaTgH8CbrT9iabaiYiYE4/sZkCDbJ3ja41osgdxOPB64FpJV5fPvd/2NxtsMyKius7lA4dI+r1T8RQdnvaWoS41FhC2/w/Fh4qI6KSunWKy3fp6S4PkSuqI6CcD3TvF1CkJiIjor+TDQAmIiOitrp1i6poERET0VgdnMXVKAiIi+qmDq7l2TQIiInqpuFAuCTFIpwJi+9alrPvFHo23s2RTezPJlt/TzkzfZS0uY7V8U3tt7bCxvb/AO2xob2nPHTaMtdLO8g1bWmkHYMmm9tqqTcdWc+2aTgVERESb0oMYLAEREf2UMYihEhAR0VOdXIupUxIQEdFfOcU0UAIiIvrJ3bvlaNckICKiv9KDGCgBERH9lXwYKAEREb2l8ZxjGiQBERH9ZHKh3BAJiIjoJeFcKDdEAiIi+isBMVACIiL6KwExUAIiIvopYxBDJSAiorcyi2mwBERE9JRzimmIBERE9JNJQAyRgIiI/soZpoESEBHRW7kOYrAERET0VwJioARERPSTDWM5xzRIAiIi+is9iIE6FRBLtoqdb2m+pOX3NN5E623tuKG9fwkt3zjWWls73LW1tbaW3bW5tba4e2MrzYzfvaGVdgDGNrf4/dWlpoCQdCbwYmCd7cdO87qAU4EXApuBN9q+spbGG7Rk1AVERIyEgXFX24Y7GzhqwOsvAB5ZbscBn51v+W1IQERETxk8Xm0bdiT7UuDOAbu8DPi8C5cBu0v6w5o+SGM6dYopIqI1ps1B6r2BX0x6vLZ87ra2CpiLBERE9Ff1MYi9JK2e9Hil7ZWzaEnTtT6L949EAiIi+qt6QKy3vWIeLa0F9p30eB/g1nkcrxUZg4iInioX66uyzd8FwBtUeCpwt+1On16ChnsQko6imNq1FDjD9oeabC8iojIDNS33LWkVcATFqai1wMnAcgDbpwHfpJji+lOKaa7H1tJwwxoLCElLgU8Dz6PoXv1Y0gW2b2iqzYiIWanpOgjbRw953cDxtTTWoiZ7EE8Gfmr7ZgBJX6KY6pWAiIgOyFIbwzQZENNN63rK1J0kHUdx4QjLdtujwXIiIiYxuMI1Dn3W5CB1pWldtlfaXmF7xbJdd22wnIiIKeq7knpRarIHsSCndUVEj2SxvoGaDIgfA4+UdADwS+A1wGsbbC8iojq7tllMi1VjAWF7u6S3ABdRTHM90/b1TbUXETFr6UEM1Oh1ELa/STH/NyKiY4zH2lu6fiHKUhsR0U8Ty33HjBIQEdFfmeY6UAIiInrJgNODGCgBERH9ZKcHMUQCIiJ6K4PUg8kdmuYl6Xbg57N8217A+gbKGbXF+LkW42eCxfm5uv6Z9rP9oPkcQNK3KT5nFettD7rn9KLUqYCYC0mr53kjj05ajJ9rMX4mWJyfazF+ppi93DAoIiKmlYCIiIhpLYaAmM2NwxeSxfi5FuNngsX5uRbjZ4pZWvBjEBER0YzF0IOIiIgGJCAiImJaCzYgJB0l6f9K+qmk9426njpI2lfSDyTdKOl6SSeMuqa6SFoq6SpJF466lrpI2l3SuZJuKv/Mnjbqmuog6e3l/3/XSVolaadR1xSjsSADQtJS4NPAC4DHAEdLesxoq6rFduCdth8NPBU4fpF8LoATgBtHXUTNTgW+bfsg4BAWweeTtDfwNmCF7cdS3MvlNaOtKkZlQQYE8GTgp7Zvtr0V+BLwshHXNG+2b7N9ZfnzRopfOHuPtqr5k7QP8CLgjFHXUhdJuwHPAv4JwPZW23eNtqraLAN2lrQM2IXcKri3FmpA7A38YtLjtSyCX6STSdofOBS4fLSV1OKTwHuAxbQy2sOB24GzylNnZ0jaddRFzZftXwIfA24BbgPutv2d0VYVo7JQA0LTPLdo5utKuh9wHnCi7Q2jrmc+JL0YWGf7ilHXUrNlwGHAZ20fCmwCFvxYmKQ9KHrjBwAPBXaVdMxoq4pRWagBsRbYd9LjfVgk3WBJyynC4Rzb54+6nhocDrxU0s8oTgU+R9IXR1tSLdYCa21P9PDOpQiMhe65wBrbt9veBpwPPH3ENcWILNSA+DHwSEkHSNqBYhDtghHXNG+SRHFO+0bbnxh1PXWwfZLtfWzvT/Hn9L9tL/h/kdr+FfALSQeWTx0J3DDCkupyC/BUSbuU/z8eySIYfI+5WZD3g7C9XdJbgIsoZlmcafv6EZdVh8OB1wPXSrq6fO79tr85wppiZm8Fzin/kXIzcOyI65k325dLOhe4kmJW3VVk2Y3eylIbERExrYV6iikiIhqWgIiIiGklICIiYloJiIiImFYCIiIippWAiM6RdLakV426joi+S0BERMS0EhAxMpL2L++jcHp5/4HvSNp5yj5HlovhXSvpTEk7ls//TNIHJV1ZvnbQaD5FxOKVgIhReyTwadsHA3cBr5x4obxRzdnAq20/juLK/7+a9N71tg8DPgu8q7WKI3oiARGjtsb2xLIiVwD7T3rtwPL1n5SPP0dxD4YJ58/wvoioQQIiRm3LpJ/H+N31waZb1n269059X0TUIAERXXYTsL+k/1A+fj1wyQjrieiVBER0lu37KFZI/bKkaynuSHfaaKuK6I+s5hoREdNKDyIiIqaVgIiIiGklICIiYloJiIiImFYCIiIippWAiIiIaSUgIiJiWv8frkwl4GoRVTYAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"dsz_zsel = dsz.isel(z_t=slice(0, k+1))\n",
"o2_zint = (dsz_zsel.O2 * dsz_zsel.dz).sum('z_t') * nmolcm2_to_molm2\n",
"o2_zint.attrs['units'] = 'mol m$^{-2}$'\n",
"o2_zint.std('time').plot()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:analysis]",
"language": "python",
"name": "conda-env-analysis-py"
},
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment