Divergence and vorticity_CvH.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 3: Divergence and vorticity\n",
    "### MAQ - 32806, Chiel van Heerwaarden & Imme Benedict, 2016\n",
    "In this assignment you will investigate whether the large-scale flow is in geostrophic and thermal wind balance. For this, you will use again the ECMWF ERA-Interim data from the 1st of January 2016 0:00. Note the new `get_pressure_index` function. This makes it easier to retrieve data at a certain pressure level.\n",
    "___\n",
    "First, we load the required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np              # Numpy is the fundamental package for scientific computing in Python.\n",
    "import netCDF4 as nc            # NetCDF is the data format of the meteorological data that we use.\n",
    "import matplotlib.pyplot as pl  # Matplotlib is a scientific plotting package.\n",
    "from mpl_toolkits.basemap import Basemap # Import the map plotting interface.\n",
    "from ipywidgets import interact\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# The statement below enforces the plots to be put into this notebook, instead of in their own windows.\n",
    "%matplotlib inline\n",
    "pl.rcParams.update({'font.size': 11})          # Set the standard font size of the plots to 11pt.\n",
    "pl.rcParams.update({'figure.figsize': [13,6]}) # Set the standard figure size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Understanding divergence, vorticity and deformation\n",
    "\n",
    "Before studying weather data, you will first study the two-dimensional flow patterns associated with divergence $\\dfrac{\\partial u}{\\partial x} + \\dfrac{\\partial v}{\\partial y}$, vorticity $\\dfrac{\\partial v}{\\partial x} - \\dfrac{\\partial u}{\\partial y}$, and the two components of deformation $\\dfrac{\\partial u}{\\partial x} - \\dfrac{\\partial v}{\\partial y}$ and $\\dfrac{\\partial v}{\\partial x} + \\dfrac{\\partial u}{\\partial y}$. Every infinitisimal velocity field can be described as a linear combination of these four. See Holton section 1.5 for a full explanation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_streamlines(div=0, zeta=0, d1=0, d2=0):\n",
    "    n = 65\n",
    "    x = np.linspace(-1., 1., n)\n",
    "    y = np.linspace(-1., 1., n)\n",
    "    u = 0.5*(div +d1)*x[np.newaxis, :] + 0.5*(d2 -zeta)*y[:,np.newaxis]\n",
    "    v = 0.5*(zeta+d2)*x[np.newaxis, :] + 0.5*(div-d1  )*y[:,np.newaxis]\n",
    "    U = (u**2 + v**2)**.5\n",
    "\n",
    "    pl.figure(figsize=(6,4))\n",
    "    pl.subplot(111, aspect='equal')\n",
    "    pl.contourf(x, y, U, 21)\n",
    "    pl.colorbar()\n",
    "    pl.streamplot(x, y, u, v, color='k', arrowsize=2.)\n",
    "    pl.xlabel('x')\n",
    "    pl.ylabel('y')\n",
    "    pl.show()\n",
    "    \n",
    "interact(plot_streamlines,\n",
    "         div=(-1,1,0.1), zeta=(-1,1,0.1), d1=(-1,1,0.1), d2=(-1,1,0.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Setting up a Python Basemap and read out 3D fields\n",
    "With the code below, you initialize the world map with the desired coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We define a map of the northern hemisphere that is stored in the m object.\n",
    "lon_start =  - 70.\n",
    "lon_end   =    30.\n",
    "lat_start =    30.\n",
    "lat_end   =    61.\n",
    "m = Basemap(llcrnrlon=lon_start,\n",
    "            llcrnrlat=lat_start,\n",
    "            urcrnrlon=lon_end,\n",
    "            urcrnrlat=lat_end,\n",
    "            projection='mill', resolution='l')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "Now, you load the data and read out the desired variables from a NetCDF file from the ECMWF ERA-Interim data archive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nc_file = nc.Dataset(\"data/era_data.nc\", \"r\")\n",
    "lat = nc_file.variables[\"latitude\"][:]\n",
    "lon = nc_file.variables[\"longitude\"][:]\n",
    "p = nc_file.variables[\"level\"][:]*100.\n",
    "t = 0\n",
    "nroll = lon.size//2\n",
    "lon = np.roll(lon, nroll)\n",
    "lon = np.where(lon>=180., lon-360., lon)\n",
    "\n",
    "u     = np.roll(nc_file.variables[\"u\"][t,:,:,:], nroll, -1)\n",
    "v     = np.roll(nc_file.variables[\"v\"][t,:,:,:], nroll, -1)\n",
    "omega = np.roll(nc_file.variables[\"w\"][t,:,:,:], nroll, -1)\n",
    "Phi   = np.roll(nc_file.variables[\"z\"][t,:,:,:], nroll, -1)\n",
    "T     = np.roll(nc_file.variables[\"t\"][t,:,:,:], nroll, -1)\n",
    "slp   = np.roll(nc_file.variables[\"msl\"][t,:,:], nroll, -1)\n",
    "pr    = np.roll(nc_file.variables[\"lsp\"][t,:,:], nroll, -1)\\\n",
    "      + np.roll(nc_file.variables[\"cp\" ][t,:,:], nroll, -1)\n",
    "\n",
    "nc_file.close()\n",
    "\n",
    "lons, lats = m(*np.meshgrid(lon, lat))\n",
    "lonrad, latrad = np.meshgrid(np.pi/180.*lon, np.pi/180.*lat)\n",
    "\n",
    "def get_lat_index(lat_plot):\n",
    "    return abs(lat-lat_plot).argmin()\n",
    "\n",
    "def get_pressure_index(p_plot):\n",
    "    return abs(p-p_plot).argmin()\n",
    "\n",
    "nq = 3\n",
    "qscale = 2000."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Plotting the temperature, geopotential and wind vectors\n",
    "As a starting point you find below the temperature and geopotential at 500 hPa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n = get_pressure_index(50000.)\n",
    "quiver = False # Set this flag to True if you want to enable the wind arrows.\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, T[n,:,:], 31, cmap=pl.cm.jet)\n",
    "m.colorbar()\n",
    "if (quiver):\n",
    "    qu = m.quiver(lons[::nq,::nq], lats[::nq,::nq],\n",
    "                  u[n,::nq,::nq], v[n,::nq, ::nq],\n",
    "                  pivot='mid', width=1.2e-3,\n",
    "                  scale=qscale)\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('T and $\\Phi$ at p = {0:.0f} hPa and SLP'.format(p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Assignment\n",
    "In this assignment you are going to analyze the structure of vorticity and divergence surrounding the cyclone at 48N, 20W. Before doing so, you will get familiar with divergence, vorticity and deformation.\n",
    "\n",
    "#### A first glimpse on vorticity, divergence and deformation\n",
    "Spatial variations in the velocity field can alter a moving air parcel in three possible ways. The variations can introduce *divergence*, *vorticity* and *deformation*. Use the applet at the top of the notebook and move the sliders to change the relative contribution of eacht of the three. The values in this applet are scaled and range between -1 and 1. The purpose of this exercise is to understand the patterns, not to understand its magnitudes. Make sure to read Holton section 1.5.\n",
    "\n",
    "1. Explain the flow pattern for a flow with only divergence. What happens if you change the sign?\n",
    "1. Explain the flow pattern for a flow with only vorticity. What happens if you change the sign?\n",
    "1. Explain the flow pattern for a flow with only deformation. What happens if you change the sign?\n",
    "1. Explain the flow pattern for a flow with divergence and vorticity. What happens if you change the relative magnitude of one of the two? Which combinations do you find in high and low pressure systems on the Northern hemisphere? Which not?\n",
    "1. Combine vorticity and one of the two deformation components and put both at the same magnitude. What do you observe? How is rotation generated in such a flow pattern?\n",
    "\n",
    "#### Vorticity on the weather map\n",
    "In this exercise, focus your answers on the surroundings of the cyclone at 48N and 20W.\n",
    "1. Calculate the relative vorticity $\\zeta \\equiv \\dfrac{\\partial v}{\\partial x} - \\dfrac{\\partial u}{\\partial y}$, and plot it at 850, 500 and 300 hPa. Explain the patterns and the differences you observe.\n",
    "2. Calculate and plot the planetary vorticity $f$ at the 500 hPa level. How does it look at 300 hPa?\n",
    "3. Plot the absolute vorticity $\\eta \\equiv \\zeta + f$ at 300 hPa. What do you observe? Which of its two components dominates?\n",
    "\n",
    "#### Divergence and vertical velocity on the weather map\n",
    "In this exercise, focus your answers on the surroundings of the cyclone at 48N and 20W.\n",
    "1. Calculate the horizontal divergence $\\dfrac{\\partial u}{\\partial x} + \\dfrac{\\partial v}{\\partial y}$ and plot it at 850, 500 and 300 hPa. Explain the patterns and the differences you observe.\n",
    "2. Plot the vertical velocity $\\omega$ at the 500 hPa level. What do you observe? How does it compare to the divergence plots?\n",
    "3. Plot the precipitation `pr` and $\\omega$ at the 500 hPa level in the same plot and discuss whether and why the two look similar.\n",
    "4. Make a vertical cross section along latitude of 48N of the temperature and plot the wind vectors in there with the `m.quiver` function. TIP: multiply $\\omega$ with -10 to be able to observe changes in the vertical velocity. Does the vertical profile corroborate your earlier explanations?\n",
    "___\n",
    "Work out the assignment below.\n",
    "____"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = 6.37e6\n",
    "dvdx = 1./(a * np.cos(latrad)) * np.gradient(v, axis=2) / np.gradient(lonrad, axis=1)\n",
    "dudy = 1./(a * np.cos(latrad)) * np.gradient(u*np.cos(latrad), axis=1) \\\n",
    "                               / np.gradient(latrad, axis=0)\n",
    "zeta = dvdx - dudy\n",
    "\n",
    "Omega = 7.292e-5\n",
    "f = 2.*Omega*np.sin(latrad)\n",
    "beta = 2.*Omega/a*np.cos(latrad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zeta_range = np.linspace(-3e-4, 3e-4, 21)\n",
    "\n",
    "# FIGURE 1\n",
    "n = get_pressure_index(85000.)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, zeta[n,:,:], zeta_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\zeta$ and $\\Phi$ at p = {0:.0f} hPa'.format(p[n] / 100));\n",
    "pl.tight_layout()\n",
    "\n",
    "# FIGURE 2\n",
    "n = get_pressure_index(50000.)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, zeta[n,:,:], zeta_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\zeta$ and $\\Phi$ at p = {0:.0f} hPa'.format(p[n] / 100));\n",
    "pl.tight_layout()\n",
    "\n",
    "# FIGURE 3\n",
    "n = get_pressure_index(30000.)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, zeta[n,:,:], zeta_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\zeta$ and $\\Phi$ at p = {0:.0f} hPa'.format(p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# FIGURE 1\n",
    "n = get_pressure_index(30000.)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, f, zeta_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('f and $\\Phi$ at p = {0:.0f} hPa'.format(p[n] / 100));\n",
    "pl.tight_layout()\n",
    "\n",
    "# FIGURE 2\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, zeta[n,:,:] + f, zeta_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\zeta$ and $\\Phi$ at p = {0:.0f} hPa'.format(p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dudx = 1./(a * np.cos(latrad)) * np.gradient(u, axis=2) / np.gradient(lonrad, axis=1)\n",
    "dvdy = 1./(a * np.cos(latrad)) * np.gradient(v*np.cos(latrad), axis=1) \\\n",
    "                               / np.gradient(latrad, axis=0)\n",
    "divh = dudx + dvdy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "div_range = np.linspace(-7e-5, 7e-5, 21)\n",
    "\n",
    "# FIGURE 1\n",
    "n = get_pressure_index(85000.)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, divh[n,:,:], div_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('div(V) at p = {0:.0f} hPa and $\\Phi$ at p = {1:.0f} hPa'.format(p[n] / 100, p[n] / 100));\n",
    "pl.tight_layout()\n",
    "\n",
    "# FIGURE 2\n",
    "n = get_pressure_index(50000.)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, divh[n,:,:], div_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('div(V) at p = {0:.0f} hPa and $\\Phi$ at p = {1:.0f} hPa'.format(p[n] / 100, p[n] / 100));\n",
    "pl.tight_layout()\n",
    "\n",
    "# FIGURE 3\n",
    "n = get_pressure_index(30000.)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, divh[n,:,:], div_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('div(V) at p = {0:.0f} hPa and $\\Phi$ at p = {1:.0f} hPa'.format(p[n] / 100, p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = get_pressure_index(50000.)\n",
    "\n",
    "omega_scale = np.linspace(-1., 1., 21)\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, omega[n,:,:], omega_scale,\n",
    "           cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 15, colors='k', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\omega$ at p = {0:.0f} hPa'.format(p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = get_pressure_index(50000.)\n",
    "\n",
    "pr_range = np.linspace(1., 8., 15)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, pr*1000, pr_range,\n",
    "           cmap=pl.cm.Greens, extend='max')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 16, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "omega_scale = np.arange(-2.1, -.2, 0.2)\n",
    "co = m.contour(lons, lats, omega[n,:,:], omega_scale, colors='r', linestyles='-', linewidths=1.5)\n",
    "\n",
    "cp = m.contour(lons, lats, slp, 16, colors='r', linestyles=':', linewidths=1.)\n",
    "pl.clabel(cp, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('dzetadt_adv_f at p = {0:.0f} hPa'.format(p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nqv = 4\n",
    "lat_plot = 48.\n",
    "j_lat = abs(lat - lat_plot).argmin()\n",
    "\n",
    "lonp, pp = np.meshgrid(lon, p)\n",
    "pl.contourf(lon, p, T[:,j_lat,:], 21)\n",
    "pl.colorbar()\n",
    "pl.quiver(lonp[:,::nqv], pp[:,::nqv], u[:,j_lat,::nqv], -10*omega[:,j_lat,::nqv], \n",
    "          pivot='mid', width=1.2e-3, scale=1000)\n",
    "pl.gca().invert_yaxis()\n",
    "pl.xlim(lon_start, lon_end)\n",
    "pl.ylim(100000, 10000);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "j_lat = get_lat_index(48.)\n",
    "nqp = 2\n",
    "qscalep = 1000.\n",
    "\n",
    "div_range = np.linspace(-7e-5, 7e-5, 29)\n",
    "pl.figure()\n",
    "\n",
    "pl.contourf(lonp, pp, divh[:,j_lat,:], div_range, cmap=pl.cm.seismic, extend='both')\n",
    "pl.colorbar()\n",
    "\n",
    "cT = pl.contour(lonp, pp, T[:,j_lat,:], 21, colors='k', linestyles='-', linewidths=1.)\n",
    "#pl.clabel(cT, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "quiver = False\n",
    "if (quiver):\n",
    "    pl.quiver(lonp[:,::nqp], pp[:,::nqp],\n",
    "              u[:,j_lat,::nqp], -10*omega[:,j_lat,::nqp],\n",
    "              pivot='mid', width=1.2e-3,\n",
    "              scale=qscalep)\n",
    "pl.gca().invert_yaxis()\n",
    "pl.xlim(lon_start, lon_end)\n",
    "pl.ylim(100000,10000)\n",
    "pl.xlabel('lon')\n",
    "pl.ylabel('p (Pa)')\n",
    "pl.title('T, u and $\\omega$ at lat = {0}'.format(lat[j_lat]))\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dzetadx = 1./(a * np.cos(latrad)) * np.gradient(zeta, axis=2) / np.gradient(lonrad, axis=1)\n",
    "dzetady = 1./a * np.gradient(zeta, axis=1) / np.gradient(latrad, axis=0)\n",
    "dzetadt_adv  = - u*dzetadx - v*dzetady\n",
    "dzetadt_beta = - v*beta\n",
    "dzetadt_div  = - (zeta + f) * divh\n",
    "dzetadt = dzetadt_adv + dzetadt_beta + dzetadt_div\n",
    "\n",
    "n500 = get_pressure_index(70000.)\n",
    "dzetadt_adv_500  = np.mean(dzetadt_adv  [0:n500,:,:], axis=0)\n",
    "dzetadt_beta_500 = np.mean(dzetadt_adv_f[0:n500,:,:], axis=0)\n",
    "dzetadt_div_500  = np.mean(dzetadt_div  [0:n500,:,:], axis=0)\n",
    "dzetadt_500      = np.mean(dzetadt      [0:n500,:,:], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "div_range = np.linspace(-2e-8, 2e-8, 21)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, dzetadt_500, div_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n500,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\partial \\zeta/\\partial t$ at p = {0:.0f} hPa and $\\Phi$ at p = {1:.0f} hPa'.format(p[n] / 100, p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# dzetadt_adv\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, dzetadt_adv_500, div_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n500,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\partial \\zeta/\\partial t$ advection at p = {0:.0f} hPa and $\\Phi$ at p = {1:.0f} hPa'.format(p[n] / 100, p[n] / 100));\n",
    "pl.tight_layout()\n",
    "\n",
    "# dzetadt_div\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, dzetadt_div_500, div_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n500,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\partial \\zeta/\\partial t$ divergence at p = {0:.0f} hPa and $\\Phi$ at p = {1:.0f} hPa'.format(p[n] / 100, p[n] / 100));\n",
    "pl.tight_layout()\n",
    "\n",
    "# dzetadt_beta\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, dzetadt_beta_500, div_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi[n500,:,:], 15, colors='k', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='#666666')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,15.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,15.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\partial \\zeta/\\partial t$ beta at p = {0:.0f} hPa and $\\Phi$ at p = {1:.0f} hPa'.format(p[n] / 100, p[n] / 100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "widgets": {
   "state": {
    "c94653e35d394d228eb24552afe33eb5": {
     "views": [
      {
       "cell_index": 3
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}