Global circulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 6: Global circulation\n",
    "### MAQ - 32806, Chiel van Heerwaarden & Imme Benedict, 2016\n",
    "In this assignment, you are going to analyze daily values of temperature, geopotential and zonal and meridional wind of the last winter (1 December 2015 to 1 March 2016) with the goal to learn the structure of the global circulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "First, we load the required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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,5]}) # Set the standard figure size."
   ]
  },
  {
   "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 =  -180.\n",
    "lon_end   =   180.\n",
    "lat_start =   -90.\n",
    "lat_end   =    90.\n",
    "m = Basemap(llcrnrlon=lon_start,\n",
    "            llcrnrlat=lat_start,\n",
    "            urcrnrlon=lon_end,\n",
    "            urcrnrlat=lat_end,\n",
    "            projection='cyl', resolution='l')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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_data3.nc\", \"r\")\n",
    "lat = nc_file.variables[\"latitude\"][:]\n",
    "lon = nc_file.variables[\"longitude\"][:]\n",
    "\n",
    "nroll = lon.size//2\n",
    "lon = np.roll(lon, nroll)\n",
    "lon = np.where(lon>=180., lon-360., lon)\n",
    "\n",
    "u_850   = np.roll(nc_file.variables[\"u\"][:,0,:,:], nroll, -1)\n",
    "v_850   = np.roll(nc_file.variables[\"v\"][:,0,:,:], nroll, -1)\n",
    "Phi_850 = np.roll(nc_file.variables[\"z\"][:,0,:,:], nroll, -1)\n",
    "T_850   = np.roll(nc_file.variables[\"t\"][:,0,:,:], nroll, -1)\n",
    "\n",
    "u_300   = np.roll(nc_file.variables[\"u\"][:,1,:,:], nroll, -1)\n",
    "v_300   = np.roll(nc_file.variables[\"v\"][:,1,:,:], nroll, -1)\n",
    "Phi_300 = np.roll(nc_file.variables[\"z\"][:,1,:,:], nroll, -1)\n",
    "T_300   = np.roll(nc_file.variables[\"t\"][:,1,:,:], nroll, -1)\n",
    "time = np.arange(0, u_850.shape[0])\n",
    "\n",
    "nc_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we define a few convenient functions, similar as in earlier assignments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lons, lats = m(*np.meshgrid(lon, lat))\n",
    "\n",
    "def get_time_index(time_plot):\n",
    "    return abs(time-time_plot).argmin()\n",
    "\n",
    "def get_lat_index(lat_plot):\n",
    "    return abs(lat-lat_plot).argmin()\n",
    "\n",
    "def get_lon_index(lon_plot):\n",
    "    if(lon_plot > 180.):\n",
    "        lon_plot = lon_plot - 360.\n",
    "    return abs(lon-lon_plot).argmin()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### World maps of temperature and geopotential\n",
    "In this assignment, you are going to analyze the structure and most important elements of the global circulation. To start, below you find two maps at 300 hPa: one with temperature and one with meridional velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Phi_range = np.arange(0, 100000, 1000)\n",
    "T_range = np.linspace(210, 245, 31)\n",
    "v_range = np.linspace(-50, 50, 31)\n",
    "\n",
    "t = 10.\n",
    "nt = get_time_index(t)\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, T_300[nt,:,:], T_range, cmap=pl.cm.inferno, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi_300[nt,:,:], Phi_range, colors='w', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='k')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,30.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,30.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('T and $\\Phi$ at p = 300 hPa, t = {0} d, index = {1}'.format(t, nt));\n",
    "pl.tight_layout()\n",
    "\n",
    "pl.figure()\n",
    "m.contourf(lons, lats, v_300[nt,:,:], v_range, cmap=pl.cm.seismic, extend='both')\n",
    "m.colorbar()\n",
    "\n",
    "cz = m.contour(lons, lats, Phi_300[nt,:,:], Phi_range, colors='w', linewidths=1.5)\n",
    "pl.clabel(cz, fontsize=10., fmt='%1.0f')\n",
    "\n",
    "m.drawcoastlines(color='k')\n",
    "# draw parallels.\n",
    "parallels = np.arange(-90.,90,30.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,30.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('v and $\\Phi$ at p = 300 hPa, t = {0} d, index = {1}'.format(t, nt));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Assignment\n",
    "#### 1. The time average of the variables\n",
    "1. Recreate the two plots, but now with the time-averaged temperature $T$, geopotential $\\Phi$ and meridional velocity $v$. Do the same at 850 hPa.\n",
    "2. Do the Rossby waves average out in time? If not, describe the location of stationary waves.\n",
    "3. How do the patterns at the Northern and Southern hemisphere compare to each other? If they are different, explain why.\n",
    "4. Create maps at 300 and 850 hPa with time mean zonal wind $u$ and contours of $\\Phi$ and explain what you see.\n",
    "___\n",
    "#### 2. The zonal and time average of the variables\n",
    "1. Calculate the time and zonal averaged of the horizontal wind components $u$ and $v$ at 850 and 300 hPa.\n",
    "2. Create plots of the averaged $u$ and $v$ with latitude at the horizontal axis and put both heights in the same plot. \n",
    "3. Explain the plot of $v$. How many circulations do you observe? What do these circulations represent? Can you see in this plot whether it is summer or winter?\n",
    "4. Explain the plot of the zonal velocity $u$. What is the most prominent feature you observe in $u$? How does this feature compare between the two heights?\n",
    "5. Calculate and plot the time and zonal average temperature at both heights. Can you link the observed pattern in $T$ to that in the zonal velocity $u$?\n",
    "___\n",
    "#### 3. Mean and eddy transport\n",
    "In this part of the assignment you will decompose the meriodional transport in that done by the mean wind and that done by baroclinic eddies. The decomposition is explained at Holton, page 329. You will do the same decomposition for the meridional transport of zonal momentum $\\overline{uv}$. This you will do step by step at the 300 hPa level. We denote time averages with an overline.\n",
    "1. Calculate the time averaged zonal and meridional wind component.\n",
    "2. Create the variables $u^\\prime$ and $v^\\prime$ that are the original data with the local temporal mean subtracted.\n",
    "3. Calculate the time averaged $\\overline{u^\\prime}$ and $\\overline{v^\\prime}$. Plot the data on a map. What do you observe?\n",
    "4. Calculate the mean transport $\\overline{u}\\,\\overline{v}$ and take its zonal average.\n",
    "5. Calculate the eddy transport $\\overline{u^\\prime v^\\prime}$ and take its zonal average.\n",
    "6. Plot the mean and eddy transport together in an xy-plot and study the Northern Hemisphere.\n",
    "7. How do the magnitudes of both terms compare? What can you conclude about the importance of baroclinic instabilities in setting the atmosphere at higher latitudes into motion?\n",
    "8. Make the same plot for temperature at 850 hPa, and explain what you observe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "Work out the assignment below.\n",
    "___"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}