diff --git a/examples/mitgcm/Test xarray netcdf.ipynb b/examples/mitgcm/Test xarray netcdf.ipynb
new file mode 100644
index 0000000..8197f03
--- /dev/null
+++ b/examples/mitgcm/Test xarray netcdf.ipynb	
@@ -0,0 +1,269 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "68667024",
+   "metadata": {},
+   "source": [
+    "# Example reading netcdf output using python xarray\n",
+    "\n",
+    "Generate file 'test.nc' with:\n",
+    "\n",
+    "    julia> include(\"MITgcm_2deg8_PO4MMcarbSCH4.jl\")\n",
+    "    julia> PALEOmodel.OutputWriters.save_netcdf(paleorun.output, \"MITgcm_PO4MMcarbSCH42deg8.nc\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d640c5a0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c9a1bd15",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import xarray as xr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a5ddcd9d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "os.getcwd()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "401e5aef",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# read top-level dataset\n",
+    "ds = xr.open_dataset(\"MITgcm_PO4MMcarbSCH42deg8.nc\")\n",
+    "ds.attrs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4b80a015",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# read ocean group (= PALEO Domain ocean) from netcdf file\n",
+    "ds_ocean = xr.open_dataset(\"MITgcm_PO4MMcarbSCH42deg8.nc\", group=\"ocean\")\n",
+    "\n",
+    "# attach z coordinates (this is not automatic!)\n",
+    "ds_ocean = ds_ocean.set_coords([\"zmid\", \"zlower\", \"zupper\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8ef1c0ef",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "ds_ocean"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11dcd1af",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# timeseries for a scalar variable\n",
+    "print(ds_ocean[\"O2_total\"])\n",
+    "ds_ocean[\"O2_total\"].plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1c315b33",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# surface P_conc at last timestep\n",
+    "P_conc = ds_ocean.P_conc\n",
+    "P_conc"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2d8967e3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# P_conc at last timestep\n",
+    "P_conc_last = ds_ocean.P_conc.isel(tmodel=6) # NB: zero-based !\n",
+    "P_conc_last"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f96a43bd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# surface P_conc at last timestep\n",
+    "P_conc_last_surface = ds_ocean.P_conc.isel(tmodel=6).isel(zt=0) # NB: zero-based !\n",
+    "P_conc_last_surface"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "937c6815",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "P_conc_last_surface.plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "245e512c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# section P_conc at last timestep\n",
+    "P_conc_last_section = ds_ocean.P_conc.isel(tmodel=6).sel(lon=200.0, method=\"nearest\") # NB: zero-based !\n",
+    "P_conc_last_section"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ecb644c7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "P_conc_last_section.plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b70a2ee3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# sections by faceting\n",
+    "\n",
+    "P_conc_last_3sections = ds_ocean.P_conc.isel(tmodel=6).sel(lon=[70.0, 200.0, 340.0], method=\"nearest\") # NB: zero-based !\n",
+    "P_conc_last_3sections"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6a4ddf13",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "P_conc_last_3sections.plot(col=\"lon\", col_wrap=3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3c910a6d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# O_conc at last timestep\n",
+    "O2_conc_last_3sections = ds_ocean.O2_conc.isel(tmodel=6).sel(lon=[70.0, 200.0, 340.0], method=\"nearest\") # NB: zero-based !\n",
+    "O2_conc_last_3sections.plot(col=\"lon\", col_wrap=3, vmin=0)\n",
+    "\n",
+    "# H2S_conc at last timestep\n",
+    "H2S_conc_last_3sections = ds_ocean.H2S_conc.isel(tmodel=6).sel(isotopelinear=\"total\").sel(lon=[70.0, 200.0, 340.0], method=\"nearest\") # NB: zero-based !\n",
+    "H2S_conc_last_3sections.plot(col=\"lon\", col_wrap=3, vmin=0)\n",
+    "\n",
+    "# H2S d34S at last timestep\n",
+    "H2S_d34S_last_3sections = ds_ocean.H2S_conc.isel(tmodel=6).sel(isotopelinear=\"delta\").sel(lon=[70.0, 200.0, 340.0], method=\"nearest\") # NB: zero-based !\n",
+    "H2S_d34S_last_3sections.plot(col=\"lon\", col_wrap=3)\n",
+    "\n",
+    "# CH4_conc at last timestep\n",
+    "CH4_conc_last_3sections = ds_ocean.CH4_conc.isel(tmodel=6).sel(isotopelinear=\"total\").sel(lon=[70.0, 200.0, 340.0], method=\"nearest\") # NB: zero-based !\n",
+    "CH4_conc_last_3sections.plot(col=\"lon\", col_wrap=3, vmin=0)\n",
+    "\n",
+    "# CH4 d13C at last timestep\n",
+    "CH4_d13C_last_3sections = ds_ocean.CH4_conc.isel(tmodel=6).sel(isotopelinear=\"delta\").sel(lon=[70.0, 200.0, 340.0], method=\"nearest\") # NB: zero-based !\n",
+    "CH4_d13C_last_3sections.plot(col=\"lon\", col_wrap=3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "80eef703",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# column at ~Pacific OMZ\n",
+    "\n",
+    "H2S_conc_PacOMZ = ds_ocean.H2S_conc.isel(tmodel=6).sel(isotopelinear=\"total\").sel(lon=260, lat=0, method=\"nearest\") # NB: zero-based !\n",
+    "\n",
+    "H2S_conc_PacOMZ"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d3f24850",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# columns from 3 lon sections at lat = 0\n",
+    "H2S_conc_depth_3lon = ds_ocean.H2S_conc.isel(tmodel=6).sel(isotopelinear=\"total\").sel(lon=[70.0, 200.0, 260.0, 340.0], lat=0, method=\"nearest\") \n",
+    "H2S_conc_depth_3lon.plot.line(y=\"zmid\", hue=\"lon\")\n",
+    "# TODO - use zlower, zupper to create a stepped plot"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "base",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}