ocean_only/flow_downslope/Understanding native output data from MOM6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This \"flow downslope\" example involves four sub-directories, `layer`, `rho`, `sigma` and `z`, in which the model is running in one of four coordinate configurations. To  use this notebook it is assumed you have run each of those experiments in place and have kept the output generated in default configuration.\n",
    "\n",
    "## CAVEAT: This is a tutorial of how to make vertical section plots which also illustrates some poor ways to plot data for comparison. Read through to the end.\n",
    "\n",
    "We will use matplotlib. The line `%pylab inline` loads all the necessary packages, including `numpy` and causes images to appear in the page. We will use `scipy`'s netcdf package to read MOM6 output. Note that this only works if MOM6 is compiled with `NETCDF=3`.\n",
    "\n",
    "To see this notebook with figures see https://gist.github.com/adcroft/dde8d3fafd77d0caaa5613e64f1d7eff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import scipy.io.netcdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accessing MOM6 output\n",
    "\n",
    "Let's open some MOM6 output files and see what variables we diagnosed. We'll use `scipy.io.netcdf_file` but you could also use the `netCDF4` package and `netCDF4.Dataset` with the same syntax. Both functions return objects with members (functions or dictionaries) such as `.variables`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "layer_file = scipy.io.netcdf_file('layer/prog.nc')\n",
    "rho_file = scipy.io.netcdf_file('rho/prog.nc')\n",
    "sigma_file = scipy.io.netcdf_file('sigma/prog.nc')\n",
    "z_file = scipy.io.netcdf_file('z/prog.nc')\n",
    "for v in layer_file.variables:\n",
    "    print(v,layer_file.variables[v].shape,layer_file.variables[v].long_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above output is the result of this `diag_table`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!head -15 layer/diag_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We actually only asked for \"e\", \"h\", \"u\", \"v\", \"salt\" and \"temp\" in the `diag_table`. The other variables (such as \"xh\", \"yq\", etc) are all 1-dimensional \"axis\" variables associated with netcdf dimensions. This is a consequence of the CF convention used in MOM6 output.\n",
    "\n",
    "### Coordinates to use for plotting\n",
    "\n",
    "Because MOM6 uses curvilinear coordinates in the horizontal, the axis variables for the horizontal dimensions are not always useful. \n",
    "\n",
    "This particular experiment *is* in Cartesian coordinates so we *can* use the variables \"xh\" for the h-point locations and \"xq\" for the u-points. In this instance we can read the 1D axis variables to use as coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use the CF dimension-variable as the horizontal coordinate\n",
    "xh = layer_file.variables['xh'][:] # This is the coordinate of the cell centers (h-points in 1D)\n",
    "xq = layer_file.variables['xq'][:] # This is the coordinate of the cell corners (u-points in 1D)\n",
    "xq = numpy.concatenate(([2*xq[0]-xq[1]],xq)) # Inserts left most edge of domain in to u-point coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to output requested in the `diag_table`, MOM6 normally writes a files with static coordinates and metrics to `ocean_geometry.nc`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "geom_file = scipy.io.netcdf_file('layer/ocean_geometry.nc')\n",
    "for v in geom_file.variables:\n",
    "    print(v,geom_file.variables[v].shape,geom_file.variables[v].long_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, you'll see some 2D coordinate variables such as `geolon` and `geolonb` that facilitates plotting when the coordinates are curvilinear. We can check that the data is Cartesian and matches the CF dimensions-variable we read earlier with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot( xh, '.', label='xh (CF 1D)');\n",
    "plt.plot( geom_file.variables['lonh'][:].T, '.', label='lonh (1D)');\n",
    "plt.plot( geom_file.variables['geolon'][:].T, '.', label='geolon (2D)');\n",
    "plt.legend(loc='lower right');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above plot just confirms that we have the same values in numerous forms of the horizontal coordinates (for this experiment). In general, to make plan view plots for a 3D configuration you will always what to use the 2D coordinate data such as that in `ocean_geometry.nc` but here we will be making vertical section plots which only require the 1D form of horizontal coordinates.\n",
    "\n",
    "### Reading bottom depth from `ocean_geometry.nc`\n",
    "\n",
    "Now let's plot the topography, which is contained in the variable `D` of `ocean_geometry.nc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot( geom_file.variables['D'][0,:]); plt.title('Depth');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Depth\" is a positive quantity for the ocean; the bottom is at height $z=-D(x,y)$. The topography visualized for this experiment is a flat shallow shelf on the left, a flat deep ocean on the right, and a linear slope in between.\n",
    "\n",
    "The \"flow downslope\" example is ostensibly a 2D configuration in x-z. One would expect the y-dimension to be equal to 1 but it is instead equal to 4. Each of the 4 j-slices of an array will contain the exact same data. This is because a dimension in MOM6 can not be reduced to below the width of a parallelization halo width, which is typically equal to 4 or more.\n",
    "\n",
    "This quickly illustrates that the model state is identical along the j-axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"mean square |h(j=0)-h(j=3)|^2 =\",\n",
    "      (( layer_file.variables['h'][-1,:,0,:]-layer_file.variables['h'][-1,:,3,:] )**2).sum() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So from here on, we will use j=0 in all plots.\n",
    "\n",
    "### Exploring vertically distributed model output\n",
    "\n",
    "Now let's look at some model data in multiple coordinate modes. We opened the output files above. The python variable `layer_file` is a handle to the netcdf file `layer/prog.nc` which is output using the traditional isopycnal (or stacked shallow water) mode of MOM6. The other models in this experiment are all ALE-mode emulating the z*-coordinate (python variable `z_file` for `z/prog.nc`), terrain-following sigma-coordinate (`sigma_file` for `sigma/prog.nc`), and continuous isopycnal coordinate (`rho_file` for `rho/prog.nc`).\n",
    "\n",
    "The diagnosed variables in each of these modes was the same. However, some axis data changes meaning. For example, the vertical coordinate in layer mode is a \"target density\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print( layer_file.variables['zl'].long_name, layer_file.variables['zl'].units, layer_file.variables['zl'][:] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the model is in ALE mode emulating a z* coordinate, then the vertical coordinate is height (although we report notional depth to aid ferret with plotting):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print( z_file.variables['zl'].long_name, z_file.variables['zl'].units, z_file.variables['zl'][:] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at salinity in the first record written by the model in each of these four coordinates. We'll plot the raw data without coordinates, i.e. in index space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "plt.subplot(221);\n",
    "plt.pcolormesh( layer_file.variables['salt'][0,:,0,:] ); plt.colorbar(); plt.title('a) Layer mode S');\n",
    "plt.subplot(222);\n",
    "plt.pcolormesh( rho_file.variables['salt'][0,:,0,:] ); plt.colorbar(); plt.title(r'b) $\\rho$-coordinate S');\n",
    "plt.subplot(223);\n",
    "plt.pcolormesh( sigma_file.variables['salt'][0,:,0,:] ); plt.colorbar(); plt.title(r'c) $\\sigma$-coordinate S');\n",
    "plt.subplot(224);\n",
    "plt.pcolormesh( z_file.variables['salt'][0,:,0,:] ); plt.colorbar(); plt.title('d) z*-coordinate S');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no topography apparent in the plots and the salinity structure is hard to make sense of! In layer mode there is not even any horizontal structure. This is because in layer mode density is homogeneous along a layer whereas in ALE mode density (in this case salinity) is allowed to vary along layers.\n",
    "\n",
    "Plotting output from any model in index-space ignores the coordinates which determines where the data is physically located. The apparent absence of topography is a symptom of this. In MOM6, layers always contain data but layers have variable thickness which can even vanish. This is what thickness looks like for the above salinity panels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "plt.subplot(221);\n",
    "plt.pcolormesh( layer_file.variables['h'][0,:,0,:] ); plt.colorbar(); plt.title('a) Layer mode h');\n",
    "plt.subplot(222);\n",
    "plt.pcolormesh( rho_file.variables['h'][0,:,0,:] ); plt.colorbar(); plt.title(r'b) $\\rho$-coordinate h');\n",
    "plt.subplot(223);\n",
    "plt.pcolormesh( sigma_file.variables['h'][0,:,0,:] ); plt.colorbar(); plt.title(r'c) $\\sigma$-coordinate h');\n",
    "plt.subplot(224);\n",
    "plt.pcolormesh( z_file.variables['h'][0,:,0,:] ); plt.colorbar(); plt.title('d) z*-coordinate h');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest distribution to explain is the terrain-following $\\sigma$-coordinate in which the layers are uniformly distributed in each column to fit the topography. Thus $h$ has no vertical structure in panel c. \n",
    "\n",
    "For the other panels it is important to remember that the k-index increase downward in the model; k=1 (fortran convention) or k=0 (python convention) is at the surface. Panel d has a region of uniform resolution (~100m) at low k which transitions to vanished thickness (~0) at some value of k in each column. You can sort of see the topography (blue region) upside down. Panels a and b look similar and have a lot of vanished regions with even surface layers being vanished.\n",
    "\n",
    "Before making sense of these thickness distributions, let's check that the total thickness in each column looks like the topography:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot( layer_file.variables['h'][0,:,0,:].sum(axis=0), label='Layer');\n",
    "plt.plot( rho_file.variables['h'][0,:,0,:].sum(axis=0), label=r'$\\rho$');\n",
    "plt.plot( sigma_file.variables['h'][0,:,0,:].sum(axis=0), label=r'$\\sigma$');\n",
    "plt.plot( z_file.variables['h'][0,:,0,:].sum(axis=0), label='z*');\n",
    "plt.legend(loc='lower right'); plt.title('Coloumn total thickness');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that although the thickness distributions are quite different between each model, the total thickness of each column reflects the topography we plotted earlier based on the `ocean_geometry.nc` file.\n",
    "\n",
    "The layer thickness almost provides enough information to calculate the actual position of quantities with which we could then make plots. The missing information is the absolute position of the top or bottom.\n",
    "\n",
    "Interfaces delineate, or bound, layers. The thickness of a layer, $h_{i,j,k}$, is related to the absolute position of the interface above, $z_{i,j,k-\\frac{1}{2}}$, and below, $z_{i,j,k+\\frac{1}{2}}$ by\n",
    "$$\n",
    "h_{i,j,k} = z_{i,j,k-\\frac{1}{2}} - z_{i,j,k+\\frac{1}{2}} \\;\\;\\; \\forall \\; k=1,2,\\ldots,nk\n",
    "$$\n",
    "\n",
    "where, by convention, integer-valued indices indicate layer-centered quantities and half-valued indices indicate interface-located quantities. Interface and and layer quantities are thus staggered in the vertical.\n",
    "\n",
    "The diagnostic variable `e` is the absolute vertical position of model interfaces. Because half-integer indices are not meaningful in most copmuter languages, there is a offset convention as follows.\n",
    "\n",
    "In FORTRAN:\n",
    "$$\n",
    "h(i,j,k) = e(i,j,k) - e(i,j,k+1) \\;\\;\\; \\forall \\; k=1,2,\\ldots,nk\n",
    "$$\n",
    "\n",
    "where arrays indices normally start at 1.\n",
    "\n",
    "In python\n",
    "$$\n",
    "h[k,j,i] = e[k,j,i] - e[k+1,j,i] \\;\\;\\; \\forall \\; k=0,1,\\ldots,nk-1\n",
    "$$\n",
    "\n",
    "where array indices start at 0. We have also indicated the [k,j,i] order of indices that arises from reading data from an model-generated netcdf file.\n",
    "\n",
    "Let's look at where the interfaces are by plotting a line for each interface (note the use of the transpose `.T` operator to get these lines plotted in the right direction):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "plt.subplot(221); plt.plot( layer_file.variables['e'][0,:,0,:].T); plt.title('a) Layer mode e');\n",
    "plt.subplot(222); plt.plot( rho_file.variables['e'][0,:,0,:].T); plt.title(r'b) $\\rho$-coordinate e');\n",
    "plt.subplot(223); plt.plot( sigma_file.variables['e'][0,:,0,:].T); plt.title(r'c) $\\sigma$-coordinate e');\n",
    "plt.subplot(224); plt.plot( z_file.variables['e'][0,:,0,:].T); plt.title('d) z*-coordinate e');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can now begin to discern the nature of the coordinates in each mode. If we zoom in on the shelf-break region we will be able to more clearly see what each coordinate mode is doing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "xl=5,12; yl=-1000,10\n",
    "plt.subplot(221); plt.plot( layer_file.variables['e'][0,:,0,:].T); plt.xlim(xl); plt.ylim(yl); plt.title('a) Layer mode e');\n",
    "plt.subplot(222); plt.plot( rho_file.variables['e'][0,:,0,:].T); plt.xlim(xl); plt.ylim(yl); plt.title(r'b) $\\rho$-coordinate e');\n",
    "plt.subplot(223); plt.plot( sigma_file.variables['e'][0,:,0,:].T); plt.xlim(xl); plt.ylim(yl); plt.title(r'c) $\\sigma$-coordinate e');\n",
    "plt.subplot(224); plt.plot( z_file.variables['e'][0,:,0,:].T); plt.xlim(xl); plt.ylim(yl); plt.title('d) z*-coordinate e');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\sigma$-coordinate (c) always has `nk` layers with finite thickness (uniformly distributed) in each column. The $z$*-coordinate model (d) seems to have a variable number of layers but in fact the layers thicknesses vanish wherever the layer would be below the topography. The isopycnal coordinates, both in layer-mode (a) and ALE-mode (b), have on one thick layer on the shelf and fewer finite-thickness layers off-shelf than the other models. In these cases, there are vanished layers at both the top and bottom of the column.\n",
    "\n",
    "So now we know the location of the interfaces we can presume the center of the layer is in between at $(e[k,j,i]+e[k+1,j,i])/2$. Let's use `contourf` to shade salinity at the layer centers. Note how we have to create a 2D \"x\" coordinate to pass to `contourf` since `contourf` expects both coordinate arrays to be 2D if either one of them is 2D. We do this by using an expression `x=xh+0*z` which uses `numpy`'s \"broadcasting\" feature (see http://docs.scipy.org/doc/numpy-1.10.1/user/basics.broadcasting.html for explanation of rules). We will also plot the interface positions on top of the shaded contours:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "xxl=50,120 # This is the zoomed-in region around the shelf break in model coordinates\n",
    "plt.subplot(221)\n",
    "z = ( layer_file.variables['e'][0,:-1,0,:] + layer_file.variables['e'][0,1:,0,:] ) / 2\n",
    "x = xh + 0*z\n",
    "plt.contourf( x, z, layer_file.variables['salt'][0,:,0,:]); plt.xlim(xxl); plt.ylim(yl); plt.title('a) Layer mode S');\n",
    "plt.plot( xh, layer_file.variables['e'][0,:,0,:].T, 'k');\n",
    "plt.subplot(222)\n",
    "z = ( rho_file.variables['e'][0,:-1,0,:] + rho_file.variables['e'][0,1:,0,:] ) / 2\n",
    "plt.contourf( x, z, rho_file.variables['salt'][0,:,0,:]); plt.xlim(xxl); plt.ylim(yl); plt.title(r'b) $\\rho$ coordinate S');\n",
    "plt.plot( xh, rho_file.variables['e'][0,:,0,:].T, 'k');\n",
    "plt.subplot(223)\n",
    "z = ( sigma_file.variables['e'][0,:-1,0,:] + sigma_file.variables['e'][0,1:,0,:] ) / 2\n",
    "plt.contourf( x, z, sigma_file.variables['salt'][0,:,0,:]); plt.xlim(xxl); plt.ylim(yl); plt.title(r'c) $\\sigma$ coordinate S');\n",
    "plt.plot( xh, sigma_file.variables['e'][0,:,0,:].T, 'k');\n",
    "plt.subplot(224)\n",
    "z = ( z_file.variables['e'][0,:-1,0,:] + z_file.variables['e'][0,1:,0,:] ) / 2\n",
    "plt.contourf( x, z, z_file.variables['salt'][0,:,0,:]); plt.xlim(xxl); plt.ylim(yl); plt.title('d) z* coordinate S');\n",
    "plt.plot( xh, z_file.variables['e'][0,:,0,:].T, 'k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above looks closer to what one imagines things look like but **there are some very big problems with the above plots**.\n",
    "\n",
    "1) The apparent topography (white regions at bottom of plots) is quite different between the panels. This happens because `contourf` only shades between cell centers and so only half of the edge cells are plotted. `contourf` does not extrapolate beyond the coordinate provided for the data location. We loose a half cell of shading at the top and bottom of the column and also at the left and right of the plot. In the isopycnal-like coordinates, the bottom layer is thick and so we loose a lot.\n",
    "\n",
    "2) The shading within layers and between columns is interpolated which is introducing interior features and gradients which should not be there. The overlaid interface positions make this apparent for the layer mode for which salinity is absolutely constant along a layer (recall first plot of salinity).\n",
    "\n",
    "To get the plot we want we essentially need to insert a layer of extra data at the top and bottom of the model and an extra column at both ends. To illustrate lets see how one might do this first for the \"layer\" output. We'll define a little function to help:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fix_contourf(nc_object, record, xh, variable='salt', clim=None, xl=None, yl=None, plot_grid=True):\n",
    "    e = nc_object.variables['e'][record,:,0,:] # Interface positions\n",
    "    z = ( e[:-1,:] + e[1:,:] ) / 2 # Layer centers\n",
    "    S = nc_object.variables[variable][record,:,0,:] # Model output\n",
    "    z = numpy.vstack( ( e[0,:], z, e[-1,:] ) ) # Add a layer at top and bottom\n",
    "    S = numpy.vstack( ( S[0,:], S, S[-1,:] ) ) # Add layer data from top and bottom\n",
    "    x = xh + 0*z\n",
    "    plt.contourf( x, z, S );\n",
    "    if clim is not None: plt.clim(clim);\n",
    "    if plot_grid: plt.plot( xh, e.T, 'k');\n",
    "    if xl is not None: plt.xlim(xl);\n",
    "    if yl is not None: plt.ylim(yl);\n",
    "\n",
    "plt.figure(figsize=(12,3))\n",
    "# Same plot as above\n",
    "plt.subplot(121)\n",
    "z = ( layer_file.variables['e'][0,:-1,0,:] + layer_file.variables['e'][0,1:,0,:] ) / 2\n",
    "x = xh + 0*z\n",
    "plt.contourf( x, z, layer_file.variables['salt'][0,:,0,:]); plt.xlim(xxl); plt.ylim(yl);\n",
    "plt.title('a) Layer mode S, as above');\n",
    "plt.plot( xh, layer_file.variables['e'][0,:,0,:].T, 'k');\n",
    "plt.clim(34,35)\n",
    "# Now with an extra layer above and below\n",
    "plt.subplot(122)\n",
    "fix_contourf(layer_file, 0, xh, xl=xxl, yl=yl, clim=(34,35)); plt.title('b) Layer mode S, plotted with extra layers');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So not the data appears to be plotted from the surface down to the topography. Using this approach for all the coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "xxl=50,120 # This is the zoomed-in region around the shelf break in model coordinates\n",
    "plt.subplot(221); fix_contourf(layer_file, 0, xh, xl=xxl, yl=yl, clim=(34,35)); plt.title('a) Layer mode S')\n",
    "plt.subplot(222); fix_contourf(rho_file, 0, xh, xl=xxl, yl=yl, clim=(34,35)); plt.title(r'b) $\\rho$ coordinate S');\n",
    "plt.subplot(223); fix_contourf(sigma_file, 0, xh, xl=xxl, yl=yl, clim=(34,35)); plt.title(r'c) $\\sigma$ coordinate S');\n",
    "plt.subplot(224); fix_contourf(z_file, 0, xh, xl=xxl, yl=yl, clim=(34,35)); plt.title('d) z* coordinate S');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A remaining issue is why does there appear to be a salinity inversion in the z*-coordinate model. Technically there is (see salinity plots in i,k-space) but the layers are vanished so we should not be seeing them. This is because `contourf` is interpolating between thick and vanished layers. The bottom-line is that `contourf` is assuming the data is smooth and interpreting data inconsistent with the model formulation which considers the data to be piecewise.\n",
    "\n",
    "### Use `pcolormesh()` to visualize\n",
    "\n",
    "The most consistent tool for visualizing piecewise data is `pcolormesh`. An important distinction between `contourf` and `pcolormesh` is that the latter takes the coordinates of the corners of cells when shading cell-centered values. Recall we loaded the coordinate `yq` and inserted an extra value on the left edge - we'll use that for the horizontal coordinate of cell edges. We will horizontally average to get an approximate position for the cell corner heights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def plot_with_pcolormesh(nc_object, record, xq, variable='salt', clim=None, xl=None, yl=None, plot_grid=True):\n",
    "    e = nc_object.variables['e'][record,:,0,:] # Interface positions for h-columns\n",
    "    ea = numpy.vstack( ( e[:,0].T, (e[:,:-1].T+e[:,1:].T)/2, e[:,-1].T ) ).T # Interface positions averaged to u-columns\n",
    "    plt.pcolormesh( xq+0*ea, ea, nc_object.variables[variable][record,:,0,:] )\n",
    "    if clim is not None: plt.clim(clim);\n",
    "    if plot_grid: plt.plot( xq, ea.T, 'k');\n",
    "    if xl is not None: plt.xlim(xl);\n",
    "    if yl is not None: plt.ylim(yl);\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "xxl=50,120 # This is the zoomed-in region around the shelf break in model coordinates\n",
    "plt.subplot(221); plot_with_pcolormesh(layer_file, 0, xq, xl=xxl, yl=yl, clim=(34,35)); plt.title('a) Layer mode S')\n",
    "plt.subplot(222); plot_with_pcolormesh(rho_file, 0, xq, xl=xxl, yl=yl, clim=(34,35)); plt.title(r'b) $\\rho$ coordinate S');\n",
    "plt.subplot(223); plot_with_pcolormesh(sigma_file, 0, xq, xl=xxl, yl=yl, clim=(34,35)); plt.title(r'c) $\\sigma$ coordinate S');\n",
    "plt.subplot(224); plot_with_pcolormesh(z_file, 0, xq, xl=xxl, yl=yl, clim=(34,35)); plt.title('d) z* coordinate S');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above plots, the vanished layers are reasonably hidden and the overall shading for salinity more similar between the plots.\n",
    "\n",
    "The above method treats each cell as a trapezoid with corners shared between neighboring cells. It does not preserve the mean depth of the cell boundaries. To give more faithful rendering of the what the model can do, a tool is provided in `MOM6-examples/tools/analysis/m6toolbox` that returns arguments that can be passed straight to `pcolormesh` consistent with various interpretations of the grid structure, e.g. pcm (piecewise constant thicknesses), plm (piecewise linear), linear (as described above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# These next two lines add the MOM6-examples/tools/analysis/ directory to the search path for python packages\n",
    "import sys\n",
    "sys.path.append('../../tools/analysis/')\n",
    "# m6toolbox is a python package that has a function that helps visualize vertical sections\n",
    "import m6toolbox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define a function to plot a section\n",
    "def plot_section(file_handle, record, xq, variable='salt', clim=None, xl=None, yl=None,  plot_grid=True, rep='pcm'):\n",
    "    \"\"\"Plots a section by reading vertical grid and scalar variable and super-sampling\n",
    "    both in order to plot vertical and horizontal reconstructions.\n",
    "    \n",
    "    Optional arguments have defaults for plotting salinity and overlaying the grid.\n",
    "    \"\"\"\n",
    "    e = file_handle.variables['e'][record,:,0,:] # Vertical grid positions\n",
    "    s = file_handle.variables[variable][record,:,0,:] # Scalar field to color\n",
    "    x,z,q = m6toolbox.section2quadmesh(xq, e, s, representation=rep) # This yields three areas at twice the model resolution\n",
    "    plt.pcolormesh(x, z, q);\n",
    "    if clim is not None: plt.clim(clim)\n",
    "    if plot_grid: plt.plot(x, z.T, 'k', hold=True);\n",
    "    if xl is not None: plt.xlim(xl)\n",
    "    if yl is not None: plt.ylim(yl)\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.subplot(2,2,1); plot_section(layer_file, 0, xq, xl=xxl, yl=yl, clim=(34,35), rep='plm'); plt.title('a) Layer S');\n",
    "plt.subplot(2,2,2); plot_section(rho_file, 0, xq, xl=xxl, yl=yl, clim=(34,35), rep='plm'); plt.title(r'b) $\\rho$-coordinate S');\n",
    "plt.subplot(2,2,3); plot_section(sigma_file, 0, xq, xl=xxl, yl=yl, clim=(34,35), rep='linear'); plt.title(r'c) $\\sigma$-coordinate S');\n",
    "plt.subplot(2,2,4); plot_section(z_file, 0, xq, xl=xxl, yl=yl, clim=(34,35), rep='pcm'); plt.title('d) z*-coordinate S');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}