{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [ENV / ATM 415: Climate Laboratory](http://www.atmos.albany.edu/facstaff/brose/classes/ENV415_Spring2016/)\n",
    "\n",
    "## Spring 2016\n",
    "\n",
    "# Lecture 13: Ice albedo feedback in the EBM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The annual mean EBM\n",
    "\n",
    "the equation is\n",
    "\n",
    "$$ C \\frac{\\partial T}{\\partial t} = (1-\\alpha) ~ Q - \\left( A + B~T \\right) + \\frac{D}{\\cos⁡\\phi } \\frac{\\partial }{\\partial \\phi} \\left(   \\cos⁡\\phi  ~ \\frac{\\partial T}{\\partial \\phi} \\right) $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Temperature-dependent ice line\n",
    "\n",
    "Let the surface albedo be larger wherever the temperature is below some threshold $T_f$:\n",
    "\n",
    "$$ \\alpha\\left(\\phi, T(\\phi) \\right) = \\left\\{\\begin{array}{ccc} \n",
    "\\alpha_0 + \\alpha_2 P_2(\\sin\\phi) & ~ & T(\\phi) > T_f \\\\\n",
    "a_i & ~ & T(\\phi) \\le T_f \\\\\n",
    "\\end{array} \\right. $$\n",
    "\n",
    "where $P_2(\\sin\\phi) = \\frac{1}{2}\\left( 3\\left(\\sin\\phi\\right)^2 - 1 \\right) $ is called the *second Legendre Polynomial* (just a mathematically convenient description of a smooth variation between the equator and pole)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import climlab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  for convenience, set up a dictionary with our reference parameters\n",
    "param = {'D':0.55, 'A':210, 'B':2, 'a0':0.3, 'a2':0.078, 'ai':0.62, 'Tf':-10.}\n",
    "model1 = climlab.EBM_annual( num_lat=180, D=0.55, A=210., B=2., Tf=-10., a0=0.3, a2=0.078, ai=0.62)\n",
    "print model1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we provided a parameter `ai` for the icy albedo, our model now contains several sub-processes contained within the process called `albedo`. Together these implement the step-function formula above.\n",
    "\n",
    "The process called `iceline` simply looks for grid cells with temperature below $T_f$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print model1.param"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  A python shortcut... we can use the dictionary to pass lots of input arguments simultaneously:\n",
    "\n",
    "#  same thing as before, but written differently:\n",
    "model1 = climlab.EBM_annual( num_lat=180, **param)\n",
    "print model1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ebm_plot(e, return_fig=False):    \n",
    "    templimits = -60,32\n",
    "    radlimits = -340, 340\n",
    "    htlimits = -6,6\n",
    "    latlimits = -90,90\n",
    "    lat_ticks = np.arange(-90,90,30)\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,12))\n",
    "\n",
    "    ax1 = fig.add_subplot(3,1,1)\n",
    "    ax1.plot(e.lat, e.Ts)\n",
    "    ax1.set_ylim(templimits)\n",
    "    ax1.set_ylabel('Temperature (deg C)')\n",
    "    \n",
    "    ax2 = fig.add_subplot(3,1,2)\n",
    "    ax2.plot(e.lat, e.ASR, 'k--', label='SW' )\n",
    "    ax2.plot(e.lat, -e.OLR, 'r--', label='LW' )\n",
    "    ax2.plot(e.lat, e.net_radiation, 'c-', label='net rad' )\n",
    "    ax2.plot(e.lat, e.heat_transport_convergence(), 'g--', label='dyn' )\n",
    "    ax2.plot(e.lat, e.net_radiation.squeeze() + e.heat_transport_convergence(), 'b-', label='total' )\n",
    "    ax2.set_ylim(radlimits)\n",
    "    ax2.set_ylabel('Energy budget (W m$^{-2}$)')\n",
    "    ax2.legend()\n",
    "    \n",
    "    ax3 = fig.add_subplot(3,1,3)\n",
    "    ax3.plot(e.lat_bounds, e.heat_transport() )\n",
    "    ax3.set_ylim(htlimits)\n",
    "    ax3.set_ylabel('Heat transport (PW)')\n",
    "    \n",
    "    for ax in [ax1, ax2, ax3]:\n",
    "        ax.set_xlabel('Latitude')\n",
    "        ax.set_xlim(latlimits)\n",
    "        ax.set_xticks(lat_ticks)\n",
    "        ax.grid()\n",
    "    \n",
    "    if return_fig:\n",
    "        return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Integrate out to equilibrium.\n",
    "model1.integrate_years(5)\n",
    "#  Check for energy balance\n",
    "print climlab.global_mean(model1.net_radiation)\n",
    "f = ebm_plot(model1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  There is a diagnostic that tells us the current location of the ice edge:\n",
    "model1.icelat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Polar-amplified warming in the EBM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add a small radiative forcing\n",
    "\n",
    "The equivalent of doubling CO2 in this model is something like \n",
    "\n",
    "$$ A \\rightarrow A - \\delta A $$\n",
    "\n",
    "where $\\delta A = 4$ W m$^{-2}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "deltaA = 4.\n",
    "\n",
    "#  This is a very handy way to \"clone\" an existing model:\n",
    "model2 = climlab.process_like(model1)\n",
    "\n",
    "#  Now change the longwave parameter:\n",
    "model2.subprocess['LW'].A = param['A'] - deltaA\n",
    "#  and integrate out to equilibrium again\n",
    "model2.integrate_years(5, verbose=False)\n",
    "\n",
    "plt.plot(model1.lat, model1.Ts)\n",
    "plt.plot(model2.lat, model2.Ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The warming is polar-amplified**:  more warming at the poles than elsewhere.\n",
    "\n",
    "Why?\n",
    "\n",
    "Also, the current ice line is now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model2.icelat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no ice left!\n",
    "\n",
    "Let's do some more greenhouse warming:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model3 = climlab.process_like(model1)\n",
    "model3.subprocess['LW'].A = param['A'] - 2*deltaA\n",
    "model3.integrate_years(5, verbose=False)\n",
    "\n",
    "plt.plot(model1.lat, model1.Ts)\n",
    "plt.plot(model2.lat, model2.Ts)\n",
    "plt.plot(model3.lat, model3.Ts)\n",
    "plt.xlim(-90, 90)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the ice-free regime, there is no polar-amplified warming. A uniform radiative forcing produces a uniform warming."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## A different kind of climate forcing: changing the solar constant\n",
    "\n",
    "Historically EBMs have been used to study the climatic response to a change in the energy output from the Sun.\n",
    "\n",
    "We can do that easily with `climlab`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = climlab.EBM_annual( num_lat=180, **param )\n",
    "#  The current (default) solar constant, corresponding to present-day conditions:\n",
    "m.subprocess.insolation.S0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  What happens if we decrease $S_0$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  First, get to equilibrium\n",
    "m.integrate_years(5.)\n",
    "#  Check for energy balance\n",
    "print climlab.global_mean(m.net_radiation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m.icelat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#  Now make the solar constant smaller:\n",
    "m.subprocess.insolation.S0 = 1300."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Integrate to new equilibrium\n",
    "m.integrate_years(10.)\n",
    "#  Check for energy balance\n",
    "print climlab.global_mean(m.net_radiation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m.icelat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ebm_plot(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "A much colder climate!  The ice line is sitting at 54º. The heat transport shows that the atmosphere is moving lots of energy across the ice line, trying hard to compensate for the strong radiative cooling everywhere poleward of the ice line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  What happens if we decrease $S_0$ even more?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Now make the solar constant smaller:\n",
    "m.subprocess.insolation.S0 = 1200.\n",
    "#  First, get to equilibrium\n",
    "m.integrate_years(5.)\n",
    "#  Check for energy balance\n",
    "print climlab.global_mean(m.net_radiation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ebm_plot(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Something very different happened! Where is the ice line now?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now what happens if we set $S_0$ back to its present-day value?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Now make the solar constant smaller:\n",
    "m.subprocess.insolation.S0 = 1365.2\n",
    "#  First, get to equilibrium\n",
    "m.integrate_years(5.)\n",
    "#  Check for energy balance\n",
    "print climlab.global_mean(m.net_radiation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ebm_plot(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is this the same climate we started with?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting out of the \"Snowball\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To melt all the ice and get out of this so-called \"Snowball Earth\", we need to raise the solar constant much higher.\n",
    "\n",
    "You will look at this in your homework assignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}