{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
" This notebook has some partially worked out example solutions -- BR "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ENV/ATM 415: Climate Laboratory\n",
"\n",
"[Brian E. J. Rose](http://www.atmos.albany.edu/facstaff/brose/index.html), University at Albany\n",
"\n",
"# Assignment 1\n",
"\n",
"## Due Thursday February 8, 2018"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Make a copy of this notebook file so you can add your answers in additional cells.\n",
"- Complete all the problems below. \n",
"- For the questions that require calculation, show your code.\n",
"- Include comments in your code to explain your method as necessary.\n",
"- Submit your solutions in a single Jupyter notebook that contains your text, your code, and your figures.\n",
"- *Try to make sure that your notebook runs cleanly without errors:*\n",
" - Save your notebook\n",
" - From the `Kernel` menu, select `Restart & Run All`\n",
" - Did the notebook run from start to finish without error and produce the expected output?\n",
" - If yes, save again and submit your notebook file\n",
" - If no, fix the errors and try again.\n",
"- Save your notebook as `[your last name].ipynb`, e.g. my notebook should be called `Rose.ipynb`. *This makes it easier for me when I collect all your answers*\n",
"- Submit your answers by email before class on **Thursday February 8 2018**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Primer section 1.6, Review question 1\n",
"\n",
"List the five reasons most persuasive, in your opinion, for building or using models. You do not need to restrict your list to climate modeling but, for each of your chosen reasons, give an example of how a model has contributed already and could help in the future.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 2: version control with `git`\n",
"\n",
"In this course we will make use of the very popular version control system called `git`. Version control means keeping track of changes to files, most often files of source code. Learning to use `git` now will simplify your workflow in the future, and help you emphasize reproducibility of your scientific computing.\n",
"\n",
"There are many resources online for learning `git`. If you have never used `git` before, take yourself through this tutorial:\n",
"\n",
"https://www.codeschool.com/courses/try-git\n",
"\n",
"It is free but will require you to register for an account on Code School.\n",
"\n",
"I also suggest reading through the *Getting Started* section here:\n",
"https://backlog.com/git-tutorial/\n",
"\n",
"You will also want to make sure that `git` is installed on your own laptop. If you are using Mac or Linux it is probably already there. (Just try typing `git` in a terminal and see what happens). Windows users may need to install if [from here](https://git-scm.com/downloads). *You will need this installed later, but it is not necessary to complete this homework assignment*.\n",
"\n",
"Once you have completed the tutorial and reading, answer these question in your own words:\n",
"\n",
"1. What was something your found surprising or especially interesting from the tutorial?\n",
"2. Why do you think it is important to use version control when writing computer code? \n",
"3. Can you think of a situation in your past experience where version control with `git` might have been helpful?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 3: Time-dependent warming in the zero-dimensional Energy Balance Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In lecture we defined a zero-dimensional energy balance model for the global mean surface temperature $T_s$ as follows\n",
"\n",
"$$ C \\frac{dT_s}{dt} = \\text{ASR} - \\text{OLR}$$\n",
"\n",
"$$ \\text{ASR} = (1-\\alpha) Q $$\n",
"\n",
"$$ \\text{OLR} = \\tau \\sigma T_s^4$$\n",
"\n",
"where we defined these terms:\n",
"\n",
"- $C$ is a heat capacity for the atmosphere-ocean column\n",
"- $\\alpha$ is the global mean planetary albedo\n",
"- $\\sigma = 5.67 \\times 10^{-8}$ W m$^{-2}$ K$^{-4}$ is the Stefan-Boltzmann constant\n",
"- $\\tau$ is our transmissivity parameter for the atmosphere.\n",
"- $Q$ is the global-mean incoming solar radiation, or *insolation*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Refer back to our class notes for parameter values.\n",
"\n",
"1. If the heat penetrated to twice as deep into the ocean, the value of $C$ would be twice as large. Would this affect the **equilibrium temperature**? Why or why not?\n",
"2. In class we used numerical timestepping to investigate a *hypothetical climate change scenario* in which $\\tau$ decreases to 0.57 and $\\alpha$ increases to 0.32. We produced a graph of $T_s(t)$ over a twenty year period, starting from an initial temperature of 288 K. Here you will repeat this calculate with a larger value of $C$ and compare the warming rates. Specifically:\n",
" - Repeat our in-class time-stepping calculation with the same parameters we used before (including a heat capacity of $C = 4\\times10^8$ J m$^{-2}$ K$^{-1}$), but extend it to 50 years. **You should create an array of temperatures with 51 elements, beginning from 288 K**.\n",
" - Now do it again, but use $C = 8\\times10^8$ J m$^{-2}$ K$^{-1}$ (representing 200 meters of water). You should **create another 51-element array** of temperatures also beginning from 288 K.\n",
" - **Make a well-labeled graph** that compares the two temperatures over the 50-year period.\n",
" \n",
"4. What do your results show about the role of heat capacity on climate change? Give a short written answer."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Closely following code we developed in Lecture 2\n",
"\n",
"def ASR(alpha,Q):\n",
" return (1-alpha)*Q\n",
"\n",
"def OLR(tau,T):\n",
" sigma = 5.67E-8\n",
" return tau*sigma*T**4\n",
"\n",
"dt = 60. * 60. * 24. * 365. # one year expressed in seconds\n",
"\n",
"c_w = 4E3 # Specific heat of water in J/kg/K\n",
"rho_w = 1E3 # Density of water in kg/m3\n",
"H = 100. # Depth of water in m\n",
"C1 = c_w * rho_w * H # Heat capacity for 100 m of water\n",
"C2 = C1 * 2. # heat capacity for 200 m of water"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Allow the user to specify a heat capacity as input argument\n",
"def step_forward(T, alpha=0.32, tau=0.57, C=C1, Q=341.3):\n",
" return T + dt / C * ( ASR(alpha,Q) - OLR(tau,T) )"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Step forward 50 years, calculate temperature, and store the result\n",
"# Do the same thing for the two different heat capacities\n",
"numsteps = 50\n",
"Tsteps_C1 = np.zeros(numsteps+1)\n",
"Tsteps_C2 = np.zeros(numsteps+1)\n",
"Years = np.zeros(numsteps+1)\n",
"Tsteps_C1[0] = 288. \n",
"Tsteps_C2[0] = 288.\n",
"for n in range(numsteps):\n",
" Years[n+1] = n+1\n",
" Tsteps_C1[n+1] = step_forward( Tsteps_C1[n], C=C1 )\n",
" Tsteps_C2[n+1] = step_forward( Tsteps_C2[n], C=C2 )"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VVX28PHvooYSEJAmQYKIDkUEQUAZFKyIDmIdy9gV\n+1j5KY6NGXXsjr6OYwEp6ggooKggWBIYFaQZBII0CQgEQXpAICHr/WOfkEtIck9yW3Lv+jzPeU65\np+xNwl05u4qqYowxxpRXlVgnwBhjTOVmgcQYY0xILJAYY4wJiQUSY4wxIbFAYowxJiQWSIwxxoTE\nAokxxpiQWCAxxhgTEgskxhhjQlIt1gmIhsMPP1xTU1PLde2uXbuoU6dOeBNUwVmeE4PlOf6Fmt95\n8+b9pqqNg52XEIEkNTWVuXPnluva9PR0+vTpE94EVXCW58RgeY5/oeZXRFb7Oc+KtowxxoTEAokx\nxpiQWCAxxhgTEgskxhhjQmKBxBhjTEgskBhjjAmJBRJjjDEhSYh+JMaU2f79sG+fW3JzD13y8tw5\neXkHb+/ff/CSn1+4LroduKgeui5Yiu4XLFDyfnHbRdZH/vwzfPdd4bHAz0vbDlTaVN1lncY7CtN+\np2ZlQVpaaDfxmc58FfbnC/u1yoHtfHVLccd8LwgasK94a+WQ/fUbNrIrtQN1UoP2KQyJBRJTee3e\nDdu2wfbth65zcgqXXbsKt3//3S179hy0fXJOjvvC3rcP9u5123HuqFgnoAQK5FKdPSQdtOyl5iHr\nfdQ4sC66vY8a5FK9yP4R5FL9wJJHtUPWRbf3U/XAduCxguPFbe+vQF+tSzJ/5g8WSExC2bcP1q+H\ndetg7Vq3rF8PmzYduvz+e/D71aoFdepA3bpuXbu2O1a/PjRr5raTkti0ZQstUlOhZk2oUaNwXaMG\nVK/ulmrVDt6uVg2qVj14u7ilShW3FGwHHhM5dFukcCm6X3SBkveL2w5YT58xg1NPPfXgz3xuq7p/\n/h07YMdOIScHdu50S8F2Tg7s2i3s2sWBJSfHxf/du931Rbf37HFLOON4zZqFPzbYR506NUr8kRas\nk6oV/lgDl4IfaXH7Jf34S/s1EDn084IfedFjRX9div7KFLfOyPiBI/t0Cd8/ZgkskJjo27kTVqyA\nZctg+fLC5eef4ddfDz2/dm1o0gQOP9ytO3SAxo3dfoMGLijUrw+HHVa4nZzsAkfVqr6StDw9nRYJ\nNHQGQF6V6mzeWYPNm2HzZtiyxa23bnUvdkXX27e7wFGwzsvz95zq1d2PomCpXbtwu3HjwthesCQl\nHbyuWdNtF10XxPviYn/BUvBlXSA9/buEGiIlP387tWtH/jkRCyQi0hIYDTQD8oE3VfVlETkeeB2o\nC2QBV6rqDhFpBHwInAiMVNU7SrhvQ2AskOpdf6mqbo1UPkyIsrNh/nz44YfCdVbWwee0aAFt28J5\n58GRR7r9lBS3tGjhAkPgt4Ep0Z49sGGD+2fPznbbGze6ZdOmg7e3bDm11OL+evVcbC6I1a1auXW9\neoXrgiU52b30JScXbhe8BLo3ARPPIvlGkgfcp6rzRSQZmCciXwDDgPtVdbqIXA8MBh4B9njrjt5S\nkgeBr1T1aRF50Nt/IIL5MH6pwpIl8PXXbvnuu4PfMI45Bnr0gEGD3HbbttCmjfu2MUHt3OlicEGJ\n3y+/FG6vXesCx7ZtxV/bqJH7679JE+jY0W3n5Kyma9dUGjVynzdsWLiuV88V1xjjR8R+VVQ1G8j2\ntneKyBKgBXAsMMM77QtgKvCIqu4CvhGRo4Pc+nygj7c9CkjHAknsrFsHU6YUBo+CwNG6NfTrByec\n4Jbjj3d/qpoS5eXB6tWulG/lSli1ygWOrCy3vWXLweeLQPPm7sWtXTs47TS337y5q/4pWDduXHxQ\nSE/Pok+f1CjkzMS7qPzNISKpQBfge2ARMAD4GLgEaFnG2zX1ghSqmi0iTcKXUuPLxo0wfjy8/z78\n73/uWLNmcMYZ7tusb18XSEyxtmyBxYshMxN++qmwimjVKteyuEBSEqSmuuXEEwu3W7Z0waN5cys2\nMhWDaITbbotIXWA68KSqThCRPwCvAI2AScBfVbVRwPnXAt1KqSPZpqqHBexvVdUGxZw3CBgE0LRp\n065jxowpV/pzcnKoW7duua6trIrLc9WcHBr/7380+fprGsyfj+Tns6tVKzaefjqbevdmd6tWlboe\nIxI/5717q/Dzz3VYsaIuq1bVISvLLVu31jhwTlLSfo444ndSUn4nJWU3LVq47RYtfqdhw30R/Se1\n3+34F2p++/btO09VuwU9UVUjtgDVcUVX95bw+THA7CLHrgVeLeWeS4Hm3nZzYGmwdHTt2lXLKy0t\nrdzXVlYH5XnNGtV771WtW9d1bzvqKNWHHlL98UfV/PyYpTHcQv0579ih+vXXqs89p3rllart26tW\nqVLYWzA5WbVHD9Xrr1d9/nnVyZNVs7JU9+8PT/rLI+F/txNAqPkF5qqP7/pIttoSYDiwRFVfDDje\nRFU3ikgV4GFcC66ymARcAzztrT8OU5JNoIUL4fnn4b//dd+Ff/4z/PWv0L17pX7zCIf9+12x1Pff\nu2XWLLdf0PchJQW6dIGLLnLrzp1dkVSC/7OZOBbJOpJewFXAQhHJ8I49BLQVkdu9/QnAiIILRCQL\nqAfUEJGBwFmqmikiw4DXVXUuLoCME5EbgDW4ehYTLt99x3EPPui+IWvXhttvh3vucW0/E1Renmu1\nnJ4O06e7aqEdO9xnDRq4hmgXXQQ9e0K3bq57izGJJJKttr4BSvob7OUSrkkt4fiNAdubgdNDTZ8p\nYtMmGDwYRo0i+bDD4Ikn4NZbXVvQBKMKixbB55+7hmjffuua3gIceyxcdhn88Y8ucBx9tL1pGGMt\nxRNdfj4MGwYPPujGrxgyhFmnnMIp/frFOmVRlZNTjQ8/dC2ZP//cjcoCrlntX/4Cp57qlmbNYptO\nYyoiCySJbMEC99Yxc6b7lvzPf6BdO/LT02OdsqhYvRomToQJE+Dbb3uRn+96cp95pusCc/bZrmO9\nMaZ0FkgSUW4u/O1v8OKLruhq9Gj3Z3cClNEsW+a6wEyYAHPnumOdOsEVV6zhllta0aOH9eg2pqzs\nv0yi2bwZLrnEzclw003wzDOuxjiOrV8P770H77zjGqOBa3z2zDNw4YWuniM9fRW9eiVugwJjQmGB\nJJEsWgQDBrhv1lGj4OqrY52iiNm9Gz7+2GXziy9cVdBJJ8HLL8MFF7je4caY8LBAkig+/tgVXyUn\nuzasPXrEOkURMX++q+oZO9a1tDrySHjoIRcz27aNdeqMiU8WSOKdKjz1FDz8sBuwaeLEuKtB3rcP\nPvwQXn3VtRuoXRsuvRSuuQZOOcVN8GOMiRwLJPEsN9d9m77/Plx5Jbz1lpspKE6sWwdvvAFvvukG\nHW7b1hVdXXONmy/DGBMdFkjiVX4+3HijCyJPPeX6icRJq6ylS+Hpp+Hdd91wJeeeC3fc4Zrt2tuH\nMdFngSRePfCAa9b797/DkCGxTk1Y/PgjPPkkfPCBG2L9ttvgrrvgqKNinTJjEpsFknj07LNuwMU7\n7nB1I5Xc7NluxJZPPnHTtz7wgBv+q4nNRGNMhWCBJN68/bb7pr3sMldhUImLs5YudVn5+GPX1eXx\nx+HOOxNy+C9jKjQLJPFk0iTXyfCss1wHikpaYbBxIwwd6irSa9WCf/zDFWHZTL3GVEwWSOLFjBlu\nzpBu3dwYIDVqBL+mgtm9G156yfU4370bbr4ZHnvMirCMqegskMSDVatcj/XUVPjsM1eRUImowrhx\ncN99rknvwIGuVdaxx8Y6ZcYYPypn2YcplJ8P117r1pMnV7pZlbKyXPPdyy6Dpk3di9XEiRZEjKlM\nSn0jEZEk4DygN3AE8DuwCPhMVRdHPnkmqH/9y337jhgBrVvHOjW+5eW5tgCPPuraA7z0kmtkZiPv\nGlP5lPjfVkQeB/4EpAPfAxuBJOAY4GkvyNynqj9GPpmmWIsXu4GkBgxw3bkriXnzXJuAH36A886D\nf//bjYlljKmcSvv7b46qPl7CZy+KSBPA/vvHSm6uG4kwOdmNEVIJmvnm5rrK82eecRXoH3zg5jqv\nBEk3xpSitECytqQPRORWVf0P7i3FxMITT7ihbidMcJULFdzKlXD55TBnDlx/PbzwgpuN0BhT+ZVW\n2T5RRLoWPSgiQ4GbIpckE9ScOW6skKuucpNrVHDvvAOdO8Py5a511vDhFkSMiSelBZJLgA9E5CQA\ncV7HVbz3iULaTHF+/90FkObN4ZVXYp2aUm3f7qZAufpq6NLFTRF/ySWxTpUxJtxKLNpS1XkiMhD3\nZnI7hW8h/VR1X1RSZw41ZIgbO+SLLyr0n/Vz5rj+kWvWuHEjH3oIqlaNdaqMMZFQ4huJiDTE1ZNc\nA7wL5AI3A3W9z0olIi1FJE1ElojIYhG5yzt+vIjMFJGFIvKJiNQLuGaIiKwQkaUicnYJ9x0pIqtE\nJMNbOpcty5XYnDmuzewdd8AZZ8Q6NSV6913o3dsN8T5jBjzyiAURY+JZaZXt8wD1tncCPYDZgHjH\ngw3enYdrHjxfRJKBeSLyBTAMuF9Vp4vI9cBg4BERaQ9cBnTA9Vn5UkSOUdX9xdx7sKp+6C+LceTh\nh6FRIze/SAW0fz/87W+uVdapp7pZCytZ/0hjTDmUVrQVUu82Vc0Gsr3tnSKyBGgBHAvM8E77ApgK\nPAKcD4xR1b3AKhFZAXQHZoaSjrjxv//BtGnw3HMVcvTCHTvcJIyffurGyHrllUo53JcxphxKK9pK\nLe1Cr/I9xc9DvHt1wXVsXAQM8D66BGjpbbcAfgm4bK13rDhPisiPIvKSiNT0k4ZKTdW9jTRr5mZz\nqmBWroSTToIpU1znwtdftyBiTCIRVS3+A5EPcIHmY1wx1yZcz/ajgb7A6cBjqvpFqQ8QqQtMB55U\n1Qki8gfgFaARMAn4q6o2EpF/AzNV9V3vuuHAZFUdX+R+zYENQA3gTWClqv69mOcOAgYBNG3atOuY\nMWN8/HMcKicnh7oxHgSxwbx5HH///Sy/807WXXhhxJ9XljxnZNTnscc6AvDYY4s54YRtkUxaxFSE\nn3O0WZ7jX6j57du37zxV7Rb0RFUtcQHaA0/ihklZCvwA/Bf4C5BU2rXe9dVxRVf3lvD5McBsb3sI\nMCTgs6nASUHu3wf4NFg6unbtquWVlpZW7mvDIj9ftWdP1ZYtVffsicoj/eZ50iTVmjVV27VTXbEi\nsmmKtJj/nGPA8hz/Qs0vMFeDfL+qaumDNqpqJvC3skSwAiIiwHBgiaq+GHC8iapuFJEqwMPA695H\nk4D/isiLuMr2trjK/aL3ba6q2d79B+KKyuLX5Mkwa5YbBqVmxSnF++9/Xf+QE05wRVqNGsU6RcaY\nWInkMPK9gKuA0wKa6vYHLheRZcBPwHpgBIC60YTHAZnA58Dt6rXYEpHJInKEd9/3RGQhsBA4HHgi\ngnmIrfx813b2qKPcUPEVxBtvuI6Gf/wjfPWVBRFjEl3EBu1W1W9wTYWL83IJ1zyJK0orerx/wPZp\nYUlgZTBxohsid9QoqF491qkB4Nln3Tzq557rBl2sVSvWKTLGxJpNbFVR7d/vhsr9wx9cu9oYU3V9\nRB54wE1CNXGiBRFjjBP0jcSri7gSOEpV/y4iRwLNVPWQ+gsTRmPHuvlGxo6NebdwVbj3XjeH1k03\nwX/+E/MkGWMqED9vJK8BJwGXe/s7gX9HLEXGTR/42GPQqRNcfHGsU8Pjj7sgctddrn7EgogxJpCf\nOpIeqnqCiPwAoKpbRcS6m0XShAmwYoUrP6oS29LHf/3LDbp4ww1uOlybhMoYU5Sfb6lcEamKN+6W\niDQG8iOaqkQ3fLibe3bAgODnRtCoUXDPPW4WwzfesCBijCmen0DyCjARaCIiTwLfABVz1MB4sHq1\nGyL+uuti+jby0UfuLeSMM+C996w4yxhTsqBFW6r6nojMww2JIsBAVV0S8ZQlqpEj3fq662KWhPnz\nD2PIEOjWzZWuVaB+kMaYCqjUQOL1Pv9RVTviOhCaSMrPhxEj3GtAq1YxScKcOfDwwx055hjXqT6B\nhiUyxpRTqWUnqpoPLPCa/JpI++orV7R1ww0xefyaNa6jYYMGuUybBg2DTl9mjDH+Wm01BxaLyGxg\nV8FBVY1tTXA8Gj7cfXsPHBj1R+/eDeefD/v2wfPP/0jz5j2ingZjTOXkJ5AMjXgqDGze7Cokbrkl\n6pUSqnD99bBgAXz2GdSq9XtUn2+Mqdz8VLZPj0ZCEt5777nXgRgUa/3zn64D/TPPwDnnQHp61JNg\njKnE/AyRspPCudtr4OYY2aWq9SKZsISi6oq1unVzvdmj6JNP3OSLV1wBgwdH9dHGmDjh543koAnC\nRWQgbi51Ey7z5sGPP7pBrKIoM9ONB9mlCwwbZh0OjTHlU+Yeb6r6EZA4Q7lHw/DhkJTkhtWNkq1b\nXeV67dqu86GN5GuMKS8/RVuBk4RXAbpRWNRlQrV7t5tu8OKL4bDDovLI/Hy4/HLX0jgtDVq2jMpj\njTFxyk+rrT8FbOcBWcD5EUlNIho/HnbsiGol+4svwtSpriStV6+oPdYYE6f8BJJhqvpt4AER6QVs\njEySEszw4dCmDZx6alQe98MP8NBDcMEFcPPNUXmkMSbO+akj+X8+j5myWrECpk93nTiiUNO9e7dr\nndW4Mbz1llWuG2PCo8Q3EhE5CTgZaCwi9wZ8VA+wsWDD4YMP3Prqq6PyuPvug59+gi+/hEaNovJI\nY0wCKK1oqwZQ1zsnsAnwDiD20/bFg8mT4YQTICUl4o+aNAlefx3uvx9OPz3ijzPGJJASA4nXo326\niIxU1dVRTFNi2LoVvvvOVVhEWHa2q8vv3BmeeCLijzPGJBg/le27ReQ5oAOQVHBQVa0vSSimTnXt\ncM89N6KPyc+Ha6+FXbtcK2ObW8QYE25+Ktvfw81F0ho3gGMWMCfYRSLSUkTSRGSJiCwWkbu848eL\nyEwRWSgin4hIvYBrhojIChFZKiJnl3Df1iLyvYgsF5GxlXb++MmTXUXFiSdG9DEvvwzTprkmv+3a\nRfRRxpgE5SeQNFLV4UCuqk5X1euBnj6uywPuU9V23vm3i0h7YBjwoKoeh5vCdzCA99lluDeffsBr\n3lzxRT0DvKSqbYGtQGwm7whFfj5MmQL9+kV0Dttly+DBB93U79bU1xgTKX4CSa63zhaRc0WkCxC0\ndlhVs1V1vre9E1gCtACOBWZ4p30BXORtnw+MUdW9qroKWEGRMb1ERHDDs3zoHRoFRH/yjlDNmQO/\n/RbRYi1VNyJ9rVrwxhvW1NcYEzl+6kieEJH6wH24/iP1gHvK8hARSQW6AN8Di4ABwMfAJUDBAB0t\ngFkBl631jgVqBGxT1bxSzqn4Jk+GKlXgrLMi9oh33nHDn7z+OjRrFrHHGGNM0DnbqwJtVfVTYDvQ\nt6wPEJG6wHjgblXdISLXA6+IyKPAJGBfwanFXF50TC8/5xQ8dxAwCKBp06akl3OSjZycnHJfW5IT\nxo5F27Xjh4ULw3rfAtu3V+fOO7vTocNu2rb9oczzi0QizxWd5TkxJFqeo5ZfVS11AdKCnVPKtdWB\nqcC9JXx+DDDb2x4CDAn4bCpwUpHzBfgNqObtnwRMDZaOrl27anmlpaWV+9pibdigCqpPPhne+wa4\n7jrVatVUf/yxfNeHPc+VgOU5MSRankPNLzBXfXzX+6kj+U5EXhWR3iJyQsES7CKvPmM4sERVXww4\n3sRbVwEeBl73PpoEXCYiNUWkNdAWmF0k6CmQRmGHyGtwRWSVx5Qpbt2/f0RuP306jBjhOh4ed1xE\nHmGMMQfxU0dysrf+e8AxJficJL2Aq4CFIpLhHXsIaCsit3v7E4ARAKq6WETGAZm4Fl+3q+p+ABGZ\nDNyoquuBB4AxIvIE8AMuWFUekydD8+Zw/PFhv/Xeva6CvXVreOSRsN/eGGOK5WeGxDLXi3jXfUPx\ndRoAL5dwzZPAk8Uc7x+w/TOVdYbG3FzXqePiiyPSjOrZZ91YWlOmuAmrjDEmGoIWbYlIUxEZLiJT\nvP32IlL5+m5UBN99B9u3R6TZ7/Ll8OST8Oc/u+4pxhgTLX7qSEbiKr6P8PaXAXdHKkFxbfJkqF49\n7KMmqsKtt7rZel96Kay3NsaYoPwEksNVdRyQD6CuD8f+iKYqXk2eDL17Q716wc8tg48+gq++gqee\nctUvxhgTTX4CyS4RaYTXX0NEeuL6lJiyWLMGFi0Ke7FWbq4bBqVdOxg0KKy3NsYYX/y02roX1zS3\njYh8CzTG5iMpu8mT3TrMzX6HDXNjak2aBNX8/DSNMSbM/LTami8ip+LGyBJgqarmBrnMFDV5smuX\ne+yxYbvlzp3w+ONwyilw3nlhu60xxpRJ0EAiIknAbcAfccVb/xOR11V1T6QTFzf27HGVGGGem/35\n52HjRvjkExuU0RgTO34KQ0YDO3EDNgJcDryDG3DR+DF9OuzeHdZirexsF0guvRS6V85eNcaYOOEn\nkByrqoHdsNNEZEGkEhSX0tJcs98+fcJ2y8cfdxXtTz0VtlsaY0y5+Gm19YPXUgsAEekBfBu5JMWh\nWbOgSxc3OUgYLFniKtlvvRXatAnLLY0xptz8BJIeuIEbs0QkC5gJnOpNlftjRFMXD/Ly3ERWPf1M\nKunPgw9C3brw8MNhu6UxxpSbn6ItG3AjFIsWufqRMAWSGTNcU9+nnoLGjcNyS2OMCYmf5r+rRaQB\nbibDagHH50cyYXFjljfpYxgCiSoMHgwtWsBdd4V8O2OMCQs/zX//AVwLrKRwNkI/w8gbcIGkSRNI\nTQ35VpMmwezZMHy4je5rjKk4/BRtXQq0UdV9Qc80h5o1y72NhNjRQ9UVZ7VuDVdfHaa0GWNMGPip\nbF8EHBbphMSlLVtg6dKwFGt99ZV7G3ngARsKxRhTsfj5SvonrgnwImBvwUFVHRCxVMWL2d5MwWEI\nJAUj+157bci3MsaYsPITSEYBzwAL8YaSNz7NmgVVqkC3biHdZuZM16fxhRegZs0wpc0YY8LETyD5\nTVVfiXhK4tGsWdCxIyQnh3SbJ5+ERo3g5pvDlC5jjAkjP3Uk80TknyJykoicULBEPGWVXX6+K9oK\nsVgrIwM++wzuvhvq1AlT2owxJoz8vJF08daB34jW/DeY5cth69aQA8k//+leaG6/PUzpMsaYMPPT\nIbFvNBISd8LQEXHpUvjgA9dSq0GDMKXLGGPCLGjRlog0FZHhIjLF228vIjdEPmmV3KxZUL9+SBNZ\nPfOMq1y/554wpssYY8LMTx3JSGAqcIS3vwy4O9hFItJSRNJEZImILBaRu7zjnUVklohkiMhcEenu\nHW8gIhNF5EcRmS0iHUu470gRWeVdnyEinf1kNOpmzXIThVTx8098qNWr4Z134KabXMd4Y4ypqPx8\nyx2uquPwmv6qah6w38d1ecB9qtoOV79yu4i0B54FhqpqZ+BRbx/gISBDVTsBVwMvl3Lvwara2Vsy\nfKQlunbtgh9/DKlY6/nnXWf4wYPDmC5jjIkAP4Fkl4g0whtny5ubZHuwi1Q1u2BgR1XdCSwBWnj3\nqeedVh9Y7223B77yzv8JSBWRpv6zUoHMnetabZUzkPz6q5tv5OqroWXLMKfNGGPCTFS19BNcU9//\nB3TEDZfSGLhEVX3PkigiqcAM7x4tcEVlggtkJ3sjDD8FJKnqvV5x13dAD1WdV+ReI4GTcL3svwIe\nVNW9FCEig4BBAE2bNu06ZswYv8k9SE5ODnXr1i3TNS3ff582b77JNx99RF79+mV+5qhRrRg5sjWj\nR39Py5a/l/n6UJUnz5Wd5TkxJFqeQ81v375956lq8B7VqlrqAtTEte7qgAsE1YGawa4LuL4uMA+4\n0Nt/BbjI274U+NLbrgeMADJwc8LPAY4v5n7NcUGoJq7X/aPB0tC1a1ctr7S0tLJfNHCgatu25Xre\nvn2qzZur9utXrsvDolx5ruQsz4kh0fIcan6Buerje95P0dZMVc1T1cWqukhVc3GzJAYlItWB8cB7\nqjrBO3wNULD9AdDdC2g7VPU6dXUnV+PefFYVvae6IjNV9xYyouD6CkO1cMTfcvjoI8jOhjvuCHO6\njDEmQkrsRyIizXDFULVEpAvuLQDcm0PQ2TBERIDhwBJVfTHgo/XAqUA6rlPjcu/8w4Dd6oarvxGY\noao7irlvc1XN9u4/EFfcVnGsWQMbNpQ7kLz6qhsqvp/NS2mMqSRK65B4Nm5CqxTgBQoDyQ5cC6tg\negFXAQtFpKBl1UPATcDLIlIN2INXjwG0A0aLyH4gEzjQV0VEJgM3qup64D0RaeylJwO4xUdaoieE\njogLF7qpdJ97DqpWDXO6jDEmQkoMJKo6ChglIhep6viy3lhVv6Ew+BTVtZjzZwJtS7hX/4Dtij00\ny6xZUKsWHHdcmS/9978hKQmuvz4C6TLGmAgJWkdSniCS0GbNcsPGV69epsu2bXMdEK+4Aho2jFDa\njDEmAsrX7doUb+9emD+/XMVao0bB7t02OKMxpvKxQBJOGRmwb1+ZA0l+Prz2Gpx0EpxgA/QbYyoZ\nX7N/i8jJQGrg+ao6OkJpqrzKWdH+5ZewbBm8+24E0mSMMREWNJCIyDtAG1wLqYIxthSwQFLUggXQ\nrBkccUTwcwP8+99uYMaLL45QuowxJoL8vJF0A9p7vRxNaRYvhg4dynRJVhZ88gn87W82H7sxpnLy\nU0eyCGgW6YRUeqqQmQnt25fpstdfdyPN23zsxpjKys8byeFApojMxg2UCICqDohYqiqjtWshJ6dM\ngeT3390ovwMHQkpKBNNmjDER5CeQPB7pRMSFzEy3LkMg+fBD2LzZmvwaYyo3P3O2T49GQiq9cgSS\nUaPgqKOgT5/IJMkYY6LBz5ztPUVkjojkiMg+EdkvIocMppjwFi+Gxo3h8MN9nf7LL/D1127yKilp\nIBljjKkE/FS2vwpcjhultxZuZN5XI5moSikzs0wttt57z9XPX3VVBNNkjDFR4Ktnu6quAKqq6n5V\nHQH0iWhQzPDRAAAckUlEQVSqKpsytthShdGjoXdvV7RljDGVmZ/K9t0iUgPIEJFngWygTmSTVclk\nZ8P27b4Dydy5sGQJvPVWhNNljDFR4OeN5CrvvDuAXUBL4KJIJqrSKWNF++jRrvPhJZdEME3GGBMl\nflptrRaRWkBzVR0ahTRVPmUIJPv2wfvvu74j9etHOF3GGBMFflpt/Qk3ztbn3n5nEZkU6YRVKosX\nu0lEmjQJeurkya7vyDXXRCFdxhgTBX6Kth4HugPbAFQ1AzcSsClQ0GLLRzve0aOhaVM488wopMsY\nY6LATyDJU9XtEU9JZaXq3kh8FGtt3gyffgpXXgnVfA3gb4wxFZ+fr7NFInIFUFVE2gJ/Bb6LbLIq\nkY0bYetWX4Fk7FjIzXWdEI0xJl74eSO5E+iAG7DxfWAHcHckE1WplKGifdQo6NQJjj8+wmkyxpgo\n8tNqazfwN28xRfkMJD/9BLNnwwsvRCFNxhgTRX5mSOwGPMShU+12ilyyKpHFi1073ubNSz3tnXfc\nvCNXXBGldBljTJT4qSN5DxgMLATy/d5YRFripuNt5l33pqq+LCKdgdeBJCAPuE1VZ4tIA+Bt3LS+\ne4DrVXVRMfdtDYwBGgLzgatUdZ/fdIWdjxZb+fkukJx9tpuJ1xhj4omfOpJNqjpJVVep6uqCxcd1\necB9qtoO6AncLiLtgWeBoaraGXjU2wf31pPhvelcDbxcwn2fAV5S1bbAVuAGH2mJHB9jbKWnu9F+\nrZLdGBOP/ASSx0RkmIhcLiIXFizBLlLVbFWd723vBJYALQAF6nmn1QfWe9vtga+8838CUkWkaeA9\nRUSA04APvUOjgIE+8hAZmza5JUggGTsW6taF88+PUrqMMSaK/BRtXQf8AahOYdGWAhP8PkREUoEu\nwPe4Fl9TReR5XCA72TttAXAh8I2IdAdaASnArwG3agRsU9U8b38tLjgV98xBwCCApk2bkp6e7je5\nB8nJySnx2voLFtAFWJCby9YSztm/Hz744GROPHEb33+fWa40RFtpeY5XlufEkGh5jlp+VbXUBVgY\n7Jwg19cF5gEXevuvABd525cCX3rb9YARuOFY3gHmAMcXuVdjYEXAfks/6evatauWV1paWskfvvaa\nKqiuWVPiKenp7pRx48qdhKgrNc9xyvKcGBItz6HmF5irPr7n/byRzBKR9qpa5j+nRaQ6MB54T1UL\n3mCuAe7ytj8AhnkBbQfu7aegCGuVtwT6DThMRKqpeytJobBoLPoyMyE5GVJSSjxl/HhISoJzzoli\nuowxJor81JH8ETcXyVIR+VFEForIj8Eu8oLBcGCJqr4Y8NF64FRv+zTczIuIyGHevCfgZmGc4QWX\nA7wImQZc7B26BvjYRx4io6CivYQWW/n5MGGCa61Vt26U02aMMVHi542kXznv3Qs3l8lCEcnwjj0E\n3AS8LCLVcM18B3mftQNGi8h+IJOA1lgiMhm4UVXXAw8AY0TkCeAHXLCKjczMUl815syBdevgn/+M\nYpqMMSbKfM1HUp4bq+o3QEmdK7oWc/5MoG0J9+ofsP0zbjTi2NqyBTZsKLXF1vjxbnDG886LYrqM\nMSbKfM3ZboqxZIlblxBIVF2x1umnQ4MGUUyXMcZEmQWS8lq82K1LCCQ//ggrV8JFNimxMSbOWSAp\nr8xMqFMHjjyy2I/Hj3dja1knRGNMvCuxjkREduI6Hh7yEa4BVb1iPkscmZnQrp2LFsUYPx569/Y1\n+64xxlRqJQYSVU2OZkIqncxMVwFSjJ9+ch+/8kqU02SMMTHge8JXEWmCG7EXAFVdE5EUVQbbt7t2\nvSXUj0zwul5ecEEU02SMMTEStI5ERAaIyHJcL/PpQBYwJcLpqtiCtNgaPx569Ci1w7sxxsQNP5Xt\n/8ANA79MVVsDpwPfRjRVFV0pLbZWrYL58621ljEmcfgJJLmquhmoIiJVVDUN6BzhdFVsmZlQqxak\nph7y0cSJbn1h0IH2jTEmPvipI9kmInWB/wHvichG3KRViWvFCmjTBqpWPeSj8ePh+OPdx8YYkwj8\nvJGcD/yOm0fkc2Al8KdIJqrCy8qC1q0PObx+PXz3nRVrGWMSS9BAoqq7cPOA9Ae2AOO8oq7EtXo1\ntGp1yOGPPnJrCyTGmETip9XWjcBs3OyFF+PmJ7k+0gmrsLZtc81/i6kf+fRTOPpo10/RGGMShZ86\nksFAl4K3EBFpBHwHvB3JhFVYWVluXSSQ/P47pKfDjTeWOD2JMcbEJT91JGuBnQH7O4FfIpOcSqCE\nQDJjhgsmNhOiMSbRlDbW1r3e5jrgexH5GDf21vm4oq7EVEIg+fxzqFkTTj31kCuMMSaulVa0VTDW\n1kpvKRC7qW0rgqwsN29uw4YHHZ4yBfr0gdq1Y5IqY4yJmdIGbRwauC8iye6w5kQ8VRVZVpZ7Gwmo\nCFm1CpYuhVtuiVmqjDEmZvy02uooIj8Ai4DFIjJPRDpEPmkVVDFNfz//3K2tfsQYk4j8VLa/Cdyr\nqq1UtRVwH/BWZJNVgRW8kQSYMsX1TzzmmJikyBhjYspPIKnjja8FgKqmA3UilqKKbNs2twQEkr17\n4euvoV8/a/ZrjElMfvqR/CwijwDvePt/wQ0pn3hWr3brgEDyzTewa5cVaxljEpefN5LrcUOkTAAm\netvXBbtIRFqKSJqILBGRxSJyl3e8s4jMEpEMEZkrIt294/VF5BMRWeCdX+wzRCRdRJZ612d4E25F\nRzFNfz//HGrUgL59o5YKY4ypUIK+kajqVuCv5bh3HnCfqs73WnzNE5EvgGeBoao6RUT6e/t9gNuB\nTFX9k4g0BpaKyHuquq+Ye1+pqnPLkabQFBNIpkxxc7PXrRv11BhjTIVQWofET3AdEIulqgNKu7Gq\nZgPZ3vZOEVkCtPDuWc87rT6wvuASIFlEBKiLGyCyYg1Xn5XlOoo0agTAL7+4Oa6uC/p+Zowx8au0\nN5Lnw/UQEUkFugDf44ajnyoiz+OK1k72TnsVmIQLLMnAn1U1v4RbjhCR/cB44AlVLTHghdXq1Qf1\nISlo9tuvX1SebowxFZJE+jvYmxRrOvCkqk4QkVeA6ao6XkQuBQap6hkicjHQC7gXaAN8ARyvqjuK\n3K+Fqq7zisvGA++q6uhinjsIGATQtGnTrmPGjClX+nNycqjrlVt1HTSIfQ0bsvDppwF49NEOLF2a\nzJgxs+KqxVZgnhOF5TkxJFqeQ81v375956lqt6AnqmqxC25MrdsD9r8HfvaWi0u6rsg9qgNTcf1Q\nCo5tpzCACbDD2/4M6B1w3tdA9yD3vxZ4NVg6unbtquWVlpZWuNOggeptt6mq6r59qsnJqoMGlfvW\nFdZBeU4QlufEkGh5DjW/wFz18V1fWqut/8MVNRWoCZyIqxi/NViA8uo6hgNLVPXFgI/WAwVDG54G\nLPe21wCne9c2BY71glbgPauJyOHednXgPFyP+8jbvh22bj1Q0f7dd7BzpxVrGWNMaXUkNVQ1cLj4\nb9TNSbJZRPx0SOwFXAUsFJEM79hDwE3AyyJSDdiDV/wE/AMYKSILcW8qD6jqbwAikqGqnXHBbKoX\nRKoCXxKtXvZF+pB8/jlUqwannx6VpxtjTIVVWiBpELijqncE7DYOdmNV/QYXEIrTtZjz1wNnlXCv\nzt56V3HXRkWRpr9TpkCvXlCvXolXGBOXcnNzWbt2LXv27Il1Usqsfv36LFmyJNbJiBq/+U1KSiIl\nJYXq1auX6zmlBZLvReQmVT3oL34RuZlEnI+k4I2kVSvWr4cFC8Crczcmoaxdu5bk5GRSU1ORStbK\nZOfOnSQnJwc/MU74ya+qsnnzZtauXUvr1q3L9ZzSAsk9wEcicgUw3zvWFVe8NLBcT6vMsrKgVi1o\n3JipI90hGxbFJKI9e/ZUyiBiiiciNGrUiE2bNpX7HqXNR7IROFlETgMKho3/TFW/LvfTKrOAeUim\nTYPmzeG442KdKGNiw4JIfAn15+lniJSvcU1xE5sXSFQhPR1OO81G+zXGGPA3aKOBA4Fk2TLYsMFN\nq2uMiY0NGzZw2WWX0aZNG9q3b0///v1ZtmwZAP369eOwww7jvPPOK/Ued999NzNmzADg1Vdf5eij\nj0ZE+O233w6co6r89a9/5eijj6ZTp07Mnz//wGejRo2ibdu2tG3bllGjRoU9j9u2beO1114L2/3u\nv/9+vv46Mu8EFkj82LEDtmyB1FTS090hCyTGxIaqcsEFF9CnTx9WrlxJZmYmTz31FL/++isAgwcP\n5p133in1Hlu2bGHWrFmccsopAPTq1Ysvv/ySVkVmP50yZQrLly9n+fLlvPnmm9x6660Hrh86dCjf\nf/89s2fPZujQoWzdujWs+SxPIFFV8vOLH1nqzjvv5OkItRDyMx+JCehDMv0jOOIIOPro2CbJmArh\n7rshIyP4eWXRuTP8618lfpyWlkb16tW55ZZbAi7pfGD79NNPJ73gL74SfPjhh/QL6E3cpUuXYs/7\n+OOPufrqqxERevbsybZt28jOziY9PZ0zzzyThg0bAnDmmWfy+eefc/nllx90fWpqKldccQVpaWnk\n5uby5ptvMmTIEFasWMHgwYMP5OG5555j3Lhx7N27lwsuuIChQ4fy4IMPsnLlSjp37syZZ57Jc889\nV+x5WVlZnHPOOfTt25eZM2fy0Ucf8dhjjzF37lxUlRtvvJF77rmHVq1asXnzZjZs2ECzZs1K/fcp\nKwskfniBRI9sRXq6exux+hFjYmPRokV07Rpad7Jvv/2Wiy++OOh569ato2XLlgf2U1JSWLduXYnH\ni9OyZUtmzpzJPffcw7XXXsu3337Lnj176NChA7fccgvTpk1j+fLlzJ49G1VlwIABzJgxg6effppF\nixaR4QXqks478sgjWbp0KSNGjOC1115j3rx5rFu3jkWLFrFz5072799/IC0nnHAC3377LRdddFF5\n/+mKZYHED68z4vL8NmRnW7GWMQeU8uZQkWVnZ9O4cdB+1QVj+h1EREo8XpwBA9yMG8cddxw5OTkk\nJyeTnJxMUlIS27ZtY9q0aUybNu3AW1FOTg7Lly/nyCOPPOg+pZ3XqlUrevbsCcBRRx3Fzz//zJ13\n3knfvn0ZOLCwt0aTJk1Yv3494WaBxI+sLEhKIn2hm4fEAokxsdOhQwc+/PDDkO5Rq1YtXz3zU1JS\n+OWXwpGi1q5dyxFHHEFKSspBxWdr166lTwlfDDVr1gSgSpUqB7YL9vPy8lBVhgwZws0333zQdVkF\no2l4SjuvTp3CUasaNGjAggULmDp1Km+99Raffvopb7/9NuD6ANWqVStovsvKKtv98FpsTZ8hNGsG\nbdvGOkHGJK7TTjuNvXv38tZbhYNuzJkzh+nTp/u+R7t27VixYkXQ8wYMGMDo0aNRVWbNmkX9+vVp\n3rw5Z599NtOmTWPr1q1s3bqVadOmcfbZZ5crP2effTZvv/02OTk5gCtO27hxI8nJyezcuTPoeUX9\n9ttv5Ofnc9FFF/Hwww8f1NJs2bJldOzYsVzpLI29kfiRlYW2SrX6EWMqABFh4sSJ3H333Tz99NMk\nJSWRmprKv7xitt69e/PTTz+Rk5NDSkoKw4cPP+RL/txzz+WNN97gxhtvBOCVV17h2WefZcOGDXTq\n1In+/fszbNgw+vfvz+TJkzn66KOpXbs2I0aMAKBhw4Y88sgjnHjiiQA8+uijByrey+qss85iyZIl\nnHTSSQDUrVuXd999lzZt2tCrVy86duzIOeecw3PPPVfseVWrVj3ofuvWreO6664jPz+f/Px8nnnm\nGcCNkbZixQq6dQs+vUiZ+RlrvrIvIc9H0qiRLrvsEQXV118v960qjUSbs0HV8lwWmZmZ4U1IFO3Y\nsePAdq9evXTr1q0xTE3kBeZ3woQJ+vDDD5d4bnE/V8IwH4kBqv7+O2zezPR9riLL6keMiQ8vvPAC\na9asiXUyoiYvL4/77rsvIve2oq0gam7YAED6r+1p2hSOOSbGCTLGhEWPHj1inYSouuSSSyJ2b3sj\nCSJpwwYUSF/W3OpHjDGmGBZIgkjasIGVtGHdpppWrGWMMcWwQBJE0q+/Mr3aGQCcemqQk40xJgFZ\nIAkiacMG0mufQ5Mm8Ic/xDo1xhhT8VggCaJm9gbS951s9SPGVCAlDSOfkZHBSSedRIcOHejUqRNj\nx44t8R6Bw8hfeeWVHHvssXTs2JHrr7+e3NxcwIaR981PG+HKvoTSjyQzubOC6muvlfsWlY71qUgM\nlbUfSX5+vvbs2VP/85//HDj2ww8/6IwZM3Tp0qW6bNkyVVVdt26dNmvW7KC+IgX9KjZv3qw9evQ4\ncPyzzz7T/Px8zc/P18suu0xf8/7Df/bZZ9qvXz/Nz8/XmTNnavfu3Q9c37p1a928ebNu2bJFW7du\nrVu2bAlrPletWqUdOnQo0zX5+fm6f//+A/uB/UiysrL0zDPPLPHaUPqRWPPf0uzaxcydboA0qx8x\n5lAxGEU+6DDyBY444giaNGnCpk2bOOywww76rOgw8v379z+w3b17d9auXQvYMPJ+WSApzerVpNOH\nxvX20K5dUqxTY4zB/zDys2fPZt++fbRp0+aQz0oaRj43N5d33nmHl19+GbBh5P2KWCARkZbAaKAZ\nkA+8qaovi0hn4HUgCcgDblPV2SJSH3gXONJL1/OqOqKY+3YFRgK1gMnAXd4rWNjpqizS6UOfbjmI\nWCAxpqiKOop8dnY2V111FaNGjaJKlUOrgksaRv62227jlFNOoXfv3oANI+9XJCvb84D7VLUd0BO4\nXUTaA88CQ1W1M/Cotw9wO5CpqscDfYAXRKRGMff9DzAIaOst/Yo5JyxWzdvCLxxJnzOqR+oRxpgy\n6tChA/PmzSvx8x07dnDuuefyxBNPHPhyLaq4YeSHDh3Kpk2bePHFFw8cK20Y+eKOF8fvMPIZGRlk\nZGSwYsUKbrjhhkPuU9p5xQ0j36dPH956660DA1NCJRxGXlWzVXW+t70TWAK0ABSo551WHygIjwok\niwvrdYEtuGB0gIg0B+qp6kzvLWQ0MJAImT7TxbFTz0uO1COMMWVU2jDy+/bt44ILLuDqq68udUiQ\nosPIDxs2jKlTp/L+++8f9AZjw8j7E5U6EhFJBboA3wN3A1NF5HlcIDvZO+1VYBIusCQDf1bVorPY\ntwDWBuyv9Y5FRPqSphxeZTPtOzaK1COMMWVU2jDy48aNY8aMGWzevJmRI0cCMHLkyEMq44sOI3/L\nLbfQqlWrA0O0X3jhhTz66KM2jLxPEqHqhcIHiNQFpgNPquoEEXkFmK6q40XkUmCQqp4hIhcDvYB7\ngTbAF8Dxqroj4F4nAv9U1TO8/d7A/6nqn4p57iBcERhNmzbtOmbMmDKnffL/bWHHDrjs9fL9glRW\nOTk51K1bN9bJiCrLs3/169fn6KOPjkCKIm///v0HvnjPOussxo0bd0iLrngSmN9PPvmEjIwMHnnk\nkWLPXbFiBdu3bz/oWN++feepavDI46eNcHkXoDowFbg34Nh2CgOYADu87c+A3gHnfQ10L3K/5sBP\nAfuXA28ES0fI85EkGMtzYqis/UhCEdivYtasWbpgwYIYpibyAvM7bty4UudfqZDzkXh1HcOBJar6\nYsBH64GCXhmnAcu97TXA6d61TYFjgZ8D76mq2cBOEenp3f9q4ONI5cEYE7969OhBp06dYp2MqLnk\nkksi9vYVyTqSXsBVwEIRKeiy9BBwE/CyiFQD9uAVPwH/AEaKyELcm8oDqvobgIhkqGvlBXArhc1/\np3iLMSaKVLXE5q6m8tEQqzgiFkhU9RtcQCjOIb2JVHU9cFYJ9+ocsD0XCH+zA2OML0lJSWzevJlG\njRpZMIkDqsrmzZtJSip/Xznr2W6MKZOUlBTWrl3Lpk2bYp2UMtuzZ09IX5iVjd/8JiUlkZKSUu7n\nWCAxxpRJ9erVad26dayTUS7p6ekHeoYngmjl14aRN8YYExILJMYYY0JigcQYY0xIIt6zvSIQkU3A\n6nJefjjwWxiTUxlYnhOD5Tn+hZrfVqp66DDJRSREIAmFiMxVP0MExBHLc2KwPMe/aOXXiraMMcaE\nxAKJMcaYkFggCe7NWCcgBizPicHyHP+ikl+rIzHGGBMSeyMxxhgTEgskpRCRfiKyVERWiMiDsU5P\nJIjI2yKyUUQWBRxrKCJfiMhyb90glmkMJxFpKSJpIrJERBaLyF3e8XjOc5KIzBaRBV6eh3rHW4vI\n916ex4pIjVinNdxEpKqI/CAin3r7cZ1nEckSkYUikiEic71jEf/dtkBSAhGpCvwbOAdoD1wuIu1j\nm6qIGAn0K3LsQeArVW0LfOXtx4s84D5VbQf0BG73fq7xnOe9wGmqejzQGegnIj2BZ4CXvDxvBW6I\nYRoj5S5gScB+IuS5r6p2Dmj2G/HfbQskJesOrFDVn1V1HzAGOD/GaQo7VZ0BbCly+HxglLc9ChgY\n1URFkKpmq+p8b3sn7kumBfGdZ1XVHG+3urcobmK5D73jcZVnABFJAc4Fhnn7QpznuQQR/922QFKy\nFsAvAftrvWOJoKk3G2XBrJRNYpyeiBCRVKAL8D1xnmeviCcD2Ah8AawEtqlqnndKPP5+/wv4PyDf\n229E/OdZgWkiMk9ECiYNjPjvtg0jX7LiZuyxJm5xQkTqAuOBu1V1R7xP0KSq+4HOInIYMBFoV9xp\n0U1V5IjIecBGVZ0nIn0KDhdzatzk2dNLVdeLSBPgCxH5KRoPtTeSkq0FWgbsp+Dmm08Ev4pIcwBv\nvTHG6QkrEamOCyLvqeoE73Bc57mAqm4D0nH1Q4d5U15D/P1+9wIGiEgWrlj6NNwbSjznuWCmWVR1\nI+4Phu5E4XfbAknJ5gBtvVYeNYDLgEkxTlO0TAKu8bavAT6OYVrCyisnHw4sUdUXAz6K5zw39t5E\nEJFawBm4uqE04GLvtLjKs6oOUdUUVU3F/d/9WlWvJI7zLCJ1RCS5YBs3dfkiovC7bR0SSyEi/XF/\nxVQF3lbVJ2OcpLATkfeBPrhRQn8FHgM+AsYBRwJrgEtUtWiFfKUkIn8E/gcspLDs/CFcPUm85rkT\nrpK1Ku6Px3Gq+ncROQr313pD4AfgL6q6N3YpjQyvaOt+VT0vnvPs5W2it1sN+K+qPikijYjw77YF\nEmOMMSGxoi1jjDEhsUBijDEmJBZIjDHGhMQCiTHGmJBYIDHGGBMSCyTGhIE434jIOQHHLhWRz2OZ\nLmOiwZr/GhMmItIR+AA3fldVIAPop6orQ7hntYCxoYypkCyQGBNGIvIssAuoA+xU1X+IyDXA7UAN\n4DvgDlXNF5E3gROAWsBYVf27d4+1wBu44f3/hRvK4yYgF1ioqn+JcraMKZUN2mhMeA0F5gP7gG7e\nW8oFwMmqmucFj8uA/wIPquoWb+ynNBH5UFUzvfvsUtVeACKSDbRS1X0FQ50YU5FYIDEmjFR1l4iM\nBXJUda+InAGcCMz1RhiuReH0BJeLyA24/4dH4CZQKwgkYwNuuxh4V0Q+xg1fY0yFYoHEmPDLp3Ac\nL8GN0/ZI4Aki0hY3e193Vd0mIu8CSQGn7ArYPhs4FTdB0cMi0tEbFt6YCsFabRkTWV8Cl4rI4QAi\n0khEjgTqATuBHd7Q3mcXd7E35XOKqn4NDAYaA7WjknJjfLI3EmMiSFUXishQ4EsRqYKrML8FmIsr\nxloE/Ax8W8ItqgH/9YYHrwI8400RbEyFYa22jDHGhMSKtowxxoTEAokxxpiQWCAxxhgTEgskxhhj\nQmKBxBhjTEgskBhjjAmJBRJjjDEhsUBijDEmJP8fBluPi9pWc8EAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(Years, Tsteps_C1, color='red', label='C1 (100 meters)')\n",
"plt.plot(Years, Tsteps_C2, color='blue', label='C2 (200 meters)')\n",
"plt.xlabel('Years')\n",
"plt.ylabel('Global mean temperature (K)');\n",
"plt.legend()\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 4: Albedo feedback in the Energy Balance Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For this exercise, we will introduce a new physical process into our model by **letting the planetary albedo depend on temperature**. The idea is that a warmer planet has less ice and snow at the surface, and thus a lower planetary albedo.\n",
"\n",
"Represent the ice-albedo feedback through the following formula:\n",
"\n",
"$$ \\alpha(T) = \\left\\{ \\begin{array}{ccc}\n",
"\\alpha_i & & T \\le T_i \\\\\n",
"\\alpha_o + (\\alpha_i-\\alpha_o) \\frac{(T-T_o)^2}{(T_i-T_o)^2} & & T_i < T < T_o \\\\\n",
"\\alpha_o & & T \\ge T_o \\end{array} \\right\\}$$\n",
"\n",
"with the following parameter values:\n",
"\n",
"- $\\alpha_o = 0.289$ is the albedo of a warm, ice-free planet\n",
"- $\\alpha_i = 0.7$ is the albedo of a very cold, completely ice-covered planet\n",
"- $T_o = 293$ K is the threshold temperature above which our model assumes the planet is ice-free\n",
"- $T_i = 260$ K is the threshold temperature below which our model assumes the planet is completely ice covered. \n",
"\n",
"For intermediate temperature, this formula gives a smooth variation in albedo with global mean temperature. It is tuned to reproduce the observed albedo $\\alpha = 0.299$ for $T = 288$ K. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. \n",
" - Define a Python function that implements the above albedo formula. *There is definitely more than one way to do it. It doesn't matter how you do it as long as it works!*\n",
" - Use your function to calculate albedos for a wide range on planetary temperature (e.g. from $T=250$ K to $T=300$ K.)\n",
" - Present your results (albedo as a function of global mean temperature, or $\\alpha(T)$) in a nicely labeled graph.\n",
" \n",
"2. Now investigate a climate change scenario with this new model:\n",
" - Suppose that the transmissivity decreases from 0.611 to 0.57 (same as before)\n",
" - Your task is to **calculate the new equilibrium temperature**. First, explain very briefly why you can't just solve for it analytically as we did when albedo was a fixed number.\n",
" - Instead, you will use numerical time-stepping to find the equilibrium temperature\n",
" - Repeat the procedure from Question 3 *(time-step forward for 50 years from an initial temperature of 288 K and make a graph of the results)*, but this time **use the function you defined above to compute the albedo for the current temperature**.\n",
" - Is the **new equilibrium temperature larger or smaller** than it was in the model with fixed albedo? Explain why very briefly."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# One way to implement the decision-making: using if / elif / else\n",
"\n",
"def intermediate_albedo(T, ao=0.289, ai=0.7, To=293., Ti=260.):\n",
" return ao + (ai-ao)*(T-To)**2/(Ti-To)**2\n",
"\n",
"def albedo(T, ao=0.289, ai=0.7, To=293., Ti=260.):\n",
" '''Given a single input temperature in Kelvin, return a single albedo'''\n",
" if T <= Ti:\n",
" return ai\n",
" elif Ti= To:\n",
" return ao\n",
" else:\n",
" print('There is something wrong with the input')"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYVOXZx/HvvRXYpS8sSFtEepGyoIJRiA0sEGPDbixY\nE43GRN/ktSRqjIktBmPs5Y0i0RhIxCAo2JUmgrSlCEiR3pa2lPv9Yw7ritvZ2bMz8/tc11ycc+Y5\nM/fD7O5vTnuOuTsiIiIASWEXICIiNYdCQURECikURESkkEJBREQKKRRERKSQQkFERAopFEREpJBC\nQURECikURESkUErYBVRUVlaW5+TkVGrd7du3k5GRUbUF1XDqc2JQnxPDofR5+vTp6929SVntYi4U\ncnJymDZtWqXWnTx5MgMHDqzagmo49TkxqM+J4VD6bGbLytNOu49ERKSQQkFERAopFEREpJBCQURE\nCikURESkUFRDwcwGm9kCM1tkZrcV8/zDZjYzeOSZ2eZo1iMiIqWL2impZpYMjAROAlYAU81srLvP\nPdDG3X9epP1PgV7RqkdERMoWzesU+gGL3H0JgJmNAoYBc0tofz5wZ7SKmbp0I/9cWMCMggXReosa\nKXXLXgaGXYSIxIxohkIL4Osi8yuAo4praGZtgLbAuyU8PwIYAZCdnc3kyZMrXMy4rwr49+ICWLyo\nwuvGKgcMp3HtdzksM3EOH+Xn51fqZySWqc+JoTr6HM1QsGKWeQlthwOvufu+4p509yeBJwFyc3O9\nMlf0DRwIpybYFZDr83dz7O8n8tHWBow8vXfY5VQbXemaGNTn6Ijm18cVQKsi8y2BVSW0HQ68EsVa\nElJWZjon5aTy5qzVzF21NexyRCQGRDMUpgLtzaytmaUR+cM/9uBGZtYRaAh8EsVaEtbgnFTq1Urh\noQmJdSxFRConaqHg7nuBG4DxwDxgtLvPMbPfmtnQIk3PB0a5e0m7luQQZKQaI447nInz1vL58k1h\nlyMiNVxUjz66+zh37+Du7dz93mDZHe4+tkibu9z9e9cwSNX5yYC2NMpI48G388IuRURquMQ5JSWB\nZaSncN3Adny4aD2fLN4QdjkiUoMpFBLERUe3IbteOg9NWID21IlISRQKCaJWajI3/LA9U5du4r28\ndWGXIyI1lEIhgZyX24qWDWvz4Nt52loQkWIpFBJIWkoSN57QntkrtzB+zpqwyxGRGkihkGDO7NWC\nw5tk8NCEBezbr60FEfkuhUKCSUlO4ucndiBvTT7/mVXSBeYikqgUCgnotO7N6dSsLg9PyGPPvv1h\nlyMiNYhCIQElJRm3nNyRpRt2MHra12WvICIJQ6GQoE7s3JTcNg15dOJCdhTsDbscEakhFAoJysy4\nbUgn1m7bzXMfLQ27HBGpIRQKCSw3pxEndm7KE5MXs2l7QdjliEgNoFBIcLee0on8gr08Pjlx7kgn\nIiVTKCS4js3qclbvlrzwyTJWbt4ZdjkiEjKFgvDzkzoA8MgEDa0tkugUCkKLBrW55Og2vD5jBXlr\ntoVdjoiESKEgAFw/6Agy0lJ44L+6badIIlMoCAANM9K4ZmA7Js5bw7SlG8MuR0RColCQQj8ZkEOT\nuun84b/zNbS2SIJSKEihOmkp3HhC5EY878xbG3Y5IhIChYJ8x3l9W9E2K4MHxs/X0NoiCUihIN+R\nmpzErad0JG9NPv/QYHkiCUehIN8zpFsz+rRpyIMT8ti+W4PliSQShYJ8j5nx69M6s27bbv72/pKw\nyxGRaqRQkGL1bt2Q03o058n3F/PNll1hlyMi1UShICW6bXAn9u+HB9/WBW0iiUKhICVq1agOlw3I\n4bUZK5izakvY5YhINVAoSKmuH3gE9Wunct+4ebqgTSQBKBSkVPXrpHLjCe35aNEGJi9YF3Y5IhJl\nCgUp04VHtSGncR3uHTePvfv2h12OiESRQkHKlJaSxG1DOrNobT6v6oI2kbimUJByOaVrNv1yGvHw\nhDy27doTdjkiEiUKBSmXAxe0rc8v4In3FoddjohEiUJByu3IVg0Y1vMwnv7gK1Zs2hF2OSISBQoF\nqZBfDu6EGfz+rflhlyIiUaBQkApp0aA21xzfjjdnreazJRvCLkdEqphCQSrs6uPacVj9Wtz977m6\n54JInFEoSIXVTkvm9lM7M3f1VkbrFFWRuKJQkEo5vUdz+uY05E/jF7Blp05RFYkXUQ0FMxtsZgvM\nbJGZ3VZCm3PNbK6ZzTGzl6NZj1QdM+POM7qycUcBj72zMOxyRKSKRC0UzCwZGAkMAboA55tZl4Pa\ntAduBwa4e1fgpmjVI1WvW4v6nJfbiuc/XsridflhlyMiVSCaWwr9gEXuvsTdC4BRwLCD2lwFjHT3\nTQDuvjaK9UgU3HJyR2qnJnPPf+aGXYqIVAGL1nDIZnY2MNjdrwzmLwaOcvcbirT5F5AHDACSgbvc\n/b/FvNYIYARAdnZ2n1GjRlWqpvz8fDIzMyu1bqyqjj6/9dUeXl1QwM/7pHNkk5Sovld56HNODOpz\nxQwaNGi6u+eW1S6av8FWzLKDEygFaA8MBFoCH5hZN3ff/J2V3J8EngTIzc31gQMHVqqgyZMnU9l1\nY1V19Ln/sfuZ8sj7jFkO1555HGkp4Z6/oM85MajP0RHN394VQKsi8y2BVcW0GePue9z9K2ABkZCQ\nGJKWksT/nt6ZJeu28+InS8MuR0QOQTRDYSrQ3szamlkaMBwYe1CbfwGDAMwsC+gALIliTRIlgzo2\nZWDHJjwycSFrt+4KuxwRqaSohYK77wVuAMYD84DR7j7HzH5rZkODZuOBDWY2F5gE3OruGjshBh04\nRbVg737uGzcv7HJEpJKielTQ3ccB4w5adkeRaQduDh4S49pmZXD18Yfz2LuLGN6vNUcf3jjskkSk\ngnRFs1Sp6wYeQYsGtblzzBz26NadIjFHoSBVqnZaMnee0YUFa7bxwsdLwy5HRCpIoSBV7qQu2QwK\nDjqv0UFnkZiiUJAqZ2bcNbQrBft00Fkk1igUJCraNM7gmuPbMWbmKj5ZrBPKRGKFQkGi5rqB7WjZ\nsDZ3jPlSB51FYoRCQaKmVmoyd53RlYVr83n+o6VhlyMi5aBQkKg6sUs2J3RqyiMT8/hmiw46i9R0\nCgWJujvP6Mre/c5v/zMn7FJEpAwKBYm61o3r8LMT2jNu9je8M29N2OWISCkUClItrvrB4XTIzuSO\nMXPYvntv2OWISAkUClIt0lKSuO/M7qzcvJNHJuaFXY6IlEChINUmN6cR5/drzbMfLeXLlVvCLkdE\niqFQkGp12+BONKyTxq/fmM2+/dG5FayIVJ5CQapV/Tqp3HFGF75YsYWXPlkadjkichCFglS7M3o0\n57gOTfjT23ms3rIz7HJEpAiFglQ7M+OeYd3Ys28/d4+dG3Y5IlKEQkFC0bpxHW48sT3/nfMNE+bq\n2gWRmkKhIKG56geH0zG7LneO+ZJ8XbsgUiMoFCQ0qclJ3Pfj7qzeuosH/js/7HJEBIWChKxPm4Zc\n1j+HFz9ZxpSvNoZdjkjCUyhI6G49pSOtGtXmV6/PYteefWGXI5LQFAoSujppKdz/4x58tX47D2sI\nDJFQKRSkRhhwRBbD+7biqfeXMGvF5rDLEUlYCgWpMf7ntM40qZvOL1+bRcFe3b5TJAwKBakx6tVK\n5d4fdWf+N9v46+TFYZcjkpAUClKjnNglm6FHHsZfJi1kwTfbwi5HJOEoFKTGufOMLtStlcovX5+l\nkVRFqplCQWqcxpnp3HlGF774ejNPf7Ak7HJEEopCQWqkoUcexsldsnnw7Tzy1mg3kkh1KXcomFma\nmXULHqnRLErEzLj3zO5k1krhltFfsGefzkYSqQ7lCgUzGwgsBEYCjwN5ZnZcFOsSoUnddO75UTdm\nr9zC45N0NpJIdSjvlsKDwMnufry7HwecAjwcvbJEIk7t3pxhPQ/jsXcX6r7OItWgvKGQ6u4LDsy4\nex6gXUhSLe4e2pVGGWncPHomu/dqbCSRaCpvKEwzs2fMbGDweAqYHs3CRA5oUCeNP5zdg7w1+Tw8\nYWHY5YjEtfKGwrXAHOBnwI3AXOCaaBUlcrBBHZtyfr9WPPn+YqYv0xDbItFSrlBw993u/pC7/9jd\nz3T3h919d7SLEynq16d14bAGtbll9BfsKNCd2kSiodRQMLPZZjarpEd1FSkCkJmewh/PPpKlG3bw\nh7d0pzaRaEgp4/nTg3+vD/59Kfj3QmBHVCoSKcUx7RrzkwE5PPfRUgZ1asrAjk3DLkkkrpS6peDu\ny9x9GTDA3X/p7rODx21ETkstlZkNNrMFZrbIzG4r5vnLzGydmc0MHldWviuSKH41uBMds+vyi3/M\nYkO+9mKKVKXyHmjOMLNjD8yYWX8go7QVzCyZyMVuQ4AuwPlm1qWYpq+6e8/g8XQ565EEVis1mUfP\n78nWXXv41euzcNegeSJVpbyhcAUw0syWmtlXRK5qvryMdfoBi9x9ibsXAKOAYZUvVeRbnZrV47bB\nnZg4by0vT1kedjkiccMq8i3LzOoF65R5aamZnQ0Mdvcrg/mLgaPc/YYibS4Dfg+sA/KAn7v718W8\n1ghgBEB2dnafUaNGlbvmovLz88nMzKzUurEqnvu8352Hpu0mb9M+7upfm8MyI99x4rnPJVGfE8Oh\n9HnQoEHT3T23zIbuXuYDyAaeAd4K5rsAV5SxzjnA00XmLwYeO6hNYyA9mL4GeLesWvr06eOVNWnS\npEqvG6vivc9rtuz0nneP91Mffd9379nn7vHf5+Koz4nhUPoMTPNy/L0v7+6j54HxwGHBfB5wUxnr\nrABaFZlvCaw6KJA2+LfXOzwF9ClnPSIANK1Xiz+c1YM5q7by4IQFZa8gIqUqbyhkuftoYD+Au+8F\nyhqEZirQ3szamlkaMBwYW7SBmTUvMjsUmFfOekQKndy1GRcc1Zon31/Cx4vXh12OSEwrbyhsN7PG\ngAOY2dFAqccVguC4gcgWxjxgtLvPMbPfmtnQoNnPzGyOmX1BZAiNyyrRBxF+c1pn2mZlcPOrX5Bf\noLORRCqrrIvXDriZyLf8dmb2EdAEOLusldx9HDDuoGV3FJm+Hbi93NWKlKBOWgp/Ht6LMx//iGe+\nTOK0kxwzC7sskZhT3rGPZgDHA/2Bq4Gu7q5hLqRG6daiPrcN6czna/fx/MdLwy5HJCaV985rtYjs\n3vkdcDdwfbBMpEa5fEAOvZomc9+4ecxasTnsckRiTnmPKbwIdAUeA/5C5JTUl0pdQyQEZsYV3dJp\nkpnODS9/ztZde8IuSSSmlDcUOrr7Fe4+KXiMADpEszCRyspMMx67oBcrN+/k9n/O1jAYIhVQ3lD4\nPDjjCAAzOwr4KDoliRy6Pm0a8YuTO/LmrNUaBkOkAko9+8jMZhM5DTUVuMTMlgfzbYjcfU2kxrr6\nuMP5ZMkG7v73XHq3bkjn5vXCLkmkxitrS+F04AxgMNCWyBlIA4Pp06JamcghSkoyHjr3SBrUTuX6\nl2ewfbfu1iZSlrJCYVsZD5EaLSsznUeH92Lp+u06viBSDmVdvDadyO6iA1cBHfiNsmD68CjVJVJl\njmnXmFtO7sgfxy+gd+sGXDagbdglidRYpYaCuxf+9phZI6A9oOsTJOZce3w7Pl++iXvenEf3lvXp\n06ZR2CWJ1EjlvXjtSuA94L/AXcG/d5S2jkhNkpRkPHhuT1o0rM11f5/Bum26jadIccp7SuqNQF9g\nmbsPAnoBGo5SYkr92qn89cI+bNm5h5++MoO9+/aHXZJIjVPeUNjl7rsAzCzd3ecDHaNXlkh0dDms\nHvf+qDufLtnIH9/W/RdEDlbeUVJXmFkD4F/ABDPbxEE3zBGJFWf1acmM5Zv423tL6NWqIYO7NQu7\nJJEao1yh4O5nBpN3mdkkoD6R4woiMemOM7rw5aqt3PqPL+iQncnhTRLrXr8iJSnv7qNC7v6eu491\n94JoFCRSHdJTknn8wt6kpiQx4qXpbNPAeSJAJUJBJF60aFCbkRf0Zun67dw0aib79+vCNhGFgiS0\nY9o15s4zuvDO/LU8OEEHnkXKe6BZJG5ddHQb5q7eyshJi+nUrB5nHHlY2CWJhEZbCpLwzIy7h3Yj\nt01Dbn3tC75cuSXskkRCo1AQAdJSkvjrRX1oWCeNES9OY32+rniWxKRQEAk0qZvOkxfnsmF7Adf9\n3wwK9uqKZ0k8CgWRIrq3rM8DZ/dgytKN3DHmSw21LQlHB5pFDjKsZwvy1mxj5KTF5GRlcM3x7cIu\nSaTaKBREinHLSR1ZtmEH9781nzaN6jCke/OwSxKpFtp9JFKMpCTjT+ccSe/WDbjp1Zl8vnxT2CWJ\nVAuFgkgJaqUm89QluTStl85VL07j6407wi5JJOoUCiKlaJyZznOX9WX33v1c/vxUtuzUGEkS3xQK\nImU4omld/nZRH75av53r/z6DPbo5j8QxhYJIOfQ/Iov7zuzOh4vW85s3dKqqxC+dfSRSTuf2bcXy\njTv4y6RFZNdL5+aTdfNBiT8KBZEKuOXkDqzdtos/v7uIJvVqcfHRbcIuSaRKKRREKsDMuO/M7mzI\nL+COMV+SlZGmaxgkruiYgkgFpSQn8ZcLetOrVQNuHDWTT5dsCLskkSqjUBCphNppyTxzaV9aN67D\nVS9MY97qrWGXJFIlFAoildQwI40XLu9HRnoKlz47RRe3SVxQKIgcghYNavPC5f3YtWcfFz3zGWu3\n7gq7JJFDolAQOUQdm9Xl+cv7sW7bbi58+jM2bi8IuySRSlMoiFSB3q0b8vSluSzfuINLnv2Mrbs0\nHIbEpqiGgpkNNrMFZrbIzG4rpd3ZZuZmlhvNekSiqX+7LP56UW/mr97G5c9NZUfB3rBLEqmwqIWC\nmSUDI4EhQBfgfDPrUky7usDPgM+iVYtIdflhp2weHd6LGcs3cfVL09m1Z1/YJYlUSDS3FPoBi9x9\nibsXAKOAYcW0+x3wAKAjdBIXTuvRnD+c1YMPFq7np698rgH0JKZEMxRaAF8XmV8RLCtkZr2AVu7+\nnyjWIVLtzsltxd1DuzJh7hp+pmCQGGLRGu3RzM4BTnH3K4P5i4F+7v7TYD4JeBe4zN2Xmtlk4Bfu\nPq2Y1xoBjADIzs7uM2rUqErVlJ+fT2ZmZqXWjVXqc7jGL93DK/ML6JOdzLVHppOSZFF5n5rU5+qi\nPlfMoEGDprt72cdt3T0qD+AYYHyR+duB24vM1wfWA0uDxy5gFZBb2uv26dPHK2vSpEmVXjdWqc/h\ne/qDJd7mV//xES9O9d179kXlPWpan6uD+lwxwDQvx9/uaO4+mgq0N7O2ZpYGDAfGFgmjLe6e5e45\n7p4DfAoM9WK2FERi2RXHtuWO07swfs4abnh5BgV7tStJaq6ohYK77wVuAMYD84DR7j7HzH5rZkOj\n9b4iNdHlx7blzjO68PbcNVyvYJAaLKpDZ7v7OGDcQcvuKKHtwGjWIhK2nwxoS5IZd46dw3V/n8HI\nC3uRnpIcdlki36ErmkWq0aX9c/jdsK5MnLeGK1+YpgvcpMZRKIhUs4uPyeGPZ/fgo0XrufiZKWzZ\nqSExpOZQKIiE4JzcVoy8oDezVmzm/Cc/ZX3+7rBLEgEUCiKhGdK9OU9f2pcl6/M592+fsGrzzrBL\nElEoiITp+A5NeOmKo1i3dTfnPPEJS9blh12SJDiFgkjI+uY04pURR7Nzzz7OfuITPl++KeySJIEp\nFERqgG4t6vP6tf3JTE/h/Kc+5Z15a8IuSRKUQkGkhmiblcHr1/anQ3ZdrnpxGi9/tjzskiQBKRRE\napAmddN55aqjOb5DE/7njdk89PaCA2OFiVQLhYJIDZORnsJTl+RyXm4r/vzuIm59bZaGxZBqE9Vh\nLkSkclKSk7j/rO40b1CLRyYuZPnGHTxxUR8aZaSFXZrEOW0piNRQZsZNJ3bg0eE9mfn1ZoaN/JC8\nNdvCLkvinEJBpIYb1rMFo68+hl179vPjxz/m3fk6M0miR6EgEgN6tmrA2BsGkJNVhytemMZT7y/R\nAWiJCoWCSIxoXr82o68+hiHdmnHvuHncMvoLdhbsC7ssiTMKBZEYUicthb+c35ufn9iBN2au5MzH\nP2LZhu1hlyVxRKEgEmOSkowbT2zPs5f1ZfWWXZz+2IfMXKv7MkjVUCiIxKhBHZvyn58eS+tGdXhk\nxm4efHsB+/brOIMcGoWCSAxr1agOr1/bnx+0SOGxdxdx2XNTWLdN92aQylMoiMS4WqnJXN4tjd//\nuDtTvtrIkEc/4P28dWGXJTFKoSASB8yM8/u1ZuwNx9IoI5VLnp3C78fN0/AYUmEKBZE40rFZXcZc\nfywXHtWav72/hLOf+Jil63V2kpSfQkEkztROS+beM7vzxEW9Wbp+O6f9+QNGTVmui92kXBQKInFq\ncLfmvHXTcfRo2YDb/jmbnzw/lW+27Aq7LKnhFAoicaxFg9r8/cqjuHtoVz5bspGTH36Pf85Yoa0G\nKZFCQSTOJSUZl/bP4a0bf0CH7LrcPPoLrn5pOmu3aqtBvk+hIJIgcrIyePXqY/j1qZ2ZnLeOEx58\nj5c+XcZ+XfAmRSgURBJIcpJx1XGHM/6m4+jRqj7/+68vOeuJj5n/zdawS5MaQqEgkoDaZmXwf1cc\nxcPnHcmyDTs4/c8fcv9b8zXqqigURBKVmXFmr5a8c/Px/Lh3C554bzEnPDiZMTNX6kB0AlMoiCS4\nhhlpPHD2kfzjmmNolJnGjaNmcvYTn/DF15vDLk1CoFAQEQD65jRizPXH8sBZPVi2YTvDRn7ELaO/\nYI3OUkooKWEXICI1R3KScW7fVgzp3oy/TFrEcx8u5c3Zq7i0fw7XHNeOhhlpYZcoUaYtBRH5nrq1\nUrl9SGcm3nw8Q7o158n3l3DcA5N47J2F5O/WDX3imUJBRErUunEdHj6vJ/+98TiOadeYByfkcfwD\nk3j6gyXsKFA4xCOFgoiUqWOzujx5SS5vXNefTs3rcs+b8xhw/7s8OnEhm3cUhF2eVCGFgoiUW6/W\nDfn7lUfz+rX96dOmIQ9PzGPA/e9y75tzdUA6TuhAs4hUWJ82DXn60r7M/2Yrf528mGc+/IrnP17K\nqd2bc2n/HHq1aoCZhV2mVIJCQUQqrVOzejw6vBe3nNSRZz/6itemr2DMzFX0aFmfS47J4fQezamV\nmhx2mVIBUd19ZGaDzWyBmS0ys9uKef4aM5ttZjPN7EMz6xLNekQkOlo3rsNdQ7vy6f+cwO+GdWVH\nwT5+8Y8vOOb373DX2Dl8uXKLrpKOEVHbUjCzZGAkcBKwAphqZmPdfW6RZi+7+xNB+6HAQ8DgaNUk\nItGVmZ7CxcfkcNHRbfh48QZembKcl6cs5/mPl9KpWV3O7tOSYT1b0KRuetilSgmiufuoH7DI3ZcA\nmNkoYBhQGAruXnRoxgxAXyVE4oCZMeCILAYckcWWHXsYO2sVr01fwT1vzuO+cfPo17YRp3ZvzuCu\nzWhar1bY5UoR0QyFFsDXReZXAEcd3MjMrgduBtKAH0axHhEJQf06qVx8dBsuProNC9ds49+zVjNu\n9mruGDOHO8fOIbdNQ07u0ozjOzahfdNMHaAOmUVrP5+ZnQOc4u5XBvMXA/3c/acltL8gaH9pMc+N\nAEYAZGdn9xk1alSlasrPzyczM7NS68Yq9TkxxGKfV+bvZ9o3e5n6zV5W5Ef+DjWqZXTLSqZ7VjKd\nGyWTmVZyQMRinw/VofR50KBB0909t6x20QyFY4C73P2UYP52AHf/fQntk4BN7l6/tNfNzc31adOm\nVaqmyZMnM3DgwEqtG6vU58QQ631euXkn7+et470F6/ho0Xq2BUNptG+aSW5OQ3LbNCI3pyGtG9Up\n3JKI9T5XxqH02czKFQrR3H00FWhvZm2BlcBw4IKiDcysvbsvDGZPAxYiIgmnRYPanN+vNef3a82e\nffuZ+fVmpny1kWlLN/KfWat5ZUpkT3S9Wil0OawenZvXI2nLHrJWbqFtVgYZ6Tq7vqpE7X/S3fea\n2Q3AeCAZeNbd55jZb4Fp7j4WuMHMTgT2AJuA7+06EpHEkpqcRN+cRvTNaQTA/v1O3tptTF+2iTmr\ntjJ31VZGTfmanXv28cyXHwKQlZlG60Z1yGmcQatGdWhSN52szHSa1E0jKzOdRhlp1ElLITlJxyvK\nEtV4dfdxwLiDlt1RZPrGaL6/iMS+pCSjU7N6dGpWr3DZvv3O6HGTqNe6M0s3bGf5hh0s27idT5ds\n4I2ZKylpr3haShJ10pKpk5pMWkoSZkZhTBiF099ZXoOc2HwvA6P8HtrmEpGYk5xkNM9MYmCP5t97\nbs++/WzIL2B9/m7W5e9m/bbdbNpRwM6C/ezYs5edBfvYUbCPgr37C8+Bd/dvz4d38Bp6dnxG6qao\nv4dCQUTiSmpyEs3q16JZ/fi7/mHy5MlRfw+NkioiIoUUCiIiUkihICIihRQKIiJSSKEgIiKFFAoi\nIlJIoSAiIoUUCiIiUihqo6RGi5mtA5ZVcvUsYH0VlhML1OfEoD4nhkPpcxt3b1JWo5gLhUNhZtPK\nM3RsPFGfE4P6nBiqo8/afSQiIoUUCiIiUijRQuHJsAsIgfqcGNTnxBD1PifUMQURESldom0piIhI\nKeIqFMyslZlNMrN5ZjbHzG4Mlt9lZivNbGbwOLXIOreb2SIzW2Bmp4RXfcWV1N/guZ8GfZpjZg8U\nWR6z/YVSP+NXi3y+S81sZpF14rXPPc3s06DP08ysX7DczOzPQZ9nmVnvcHtQcaX0+Ugz+8TMZpvZ\nv82sXpF1Yv1zrmVmU8zsi6DPdwfL25rZZ2a2MPg5TwuWpwfzi4Lnc6qkEHePmwfQHOgdTNcF8oAu\nwF3AL4pp3wX4AkgH2gKLgeSw+1EF/R0ETATSg+eaxkN/S+vzQW0eBO6I9z4DbwNDguWnApOLTL9F\n5O6SRwOfhd2HKuzzVOD4YPnlwO/i6HM2IDOYTgU+Cz6/0cDwYPkTwLXB9HXAE8H0cODVqqgjrrYU\n3H21u88IprcB84AWpawyDBjl7rvd/StgEdAv+pVWjVL6ey1wv7vvDp5bG6wS0/2Fsj9jMzPgXOCV\nYFE899kYvA6kAAAGiElEQVSBA9+U6wOrgulhwIse8SnQwMy+f9/KGqyUPncE3g+aTQDOCqbj4XN2\nd88PZlODhwM/BF4Llr8A/CiYHhbMEzx/QvDzf0jiKhSKCjalehFJW4Abgk3pZ82sYbCsBfB1kdVW\nUHqI1FgH9bcD8INgk/I9M+sbNIub/kKxnzHAD4A17r4wmI/nPt8E/NHMvgb+BNweNIvnPn8JDA2e\nOgdoFUzHRZ/NLDnY9bmWSOgtBja7+96gSdF+FfY5eH4L0PhQa4jLUDCzTOB14CZ33wr8FWgH9ARW\nE9m9AJHNtYPF3OlYxfQ3BWhIZNPzVmB08A0iLvoLxfb5gPP5disB4rvP1wI/d/dWwM+BZw40LWb1\neOnz5cD1ZjadyG6lggNNi1k95vrs7vvcvSfQksiWTufimgX/RqXPcRcKZpZK5Ifo7+7+TwB3XxP8\nZ+8HnuLbzcoVfPtNAyIfxCpiSHH9JdKvfwabo1OA/UTGTIn5/kKJfcbMUoAfA68WaR7Pfb4UODD9\nD+Lo5xpK/F2e7+4nu3sfIuG/OGgeF30+wN03A5OJfLFrEPxsw3f7Vdjn4Pn6wMZDfe+4CoXg2/Az\nwDx3f6jI8qL7U88ksgkKMBYYHhzFbwu0B6ZUV72HqqT+Av8ish8SM+sApBEZRCum+wul9hngRGC+\nu68osiye+7wKOD6Y/iFwYJfZWOCS4Cyko4Et7r662gquAqX8LjcN/k0CfkPkwCvEx+fcxMwaBNO1\nifw8zwMmAWcHzS4FxgTTY4N5guff9eCo8yEJ+4h7VT6AY4lsPs0CZgaPU4GXgNnB8rFA8yLr/JrI\nt40FBGdyxMqjlP6mAf9HJPxmAD+Mh/6W1ufgueeBa4pZJy77HCyfTuSsm8+APkF7A0YGfZ4N5Ibd\nhyrs841EzkTKA+4nuAA3Tj7nHsDnQZ+/5Nsz6A4nEnCLiGwRHjirsFYwvyh4/vCqqENXNIuISKG4\n2n0kIiKHRqEgIiKFFAoiIlJIoSAiIoUUCiIiUkihIDHBzBoXGQX1G/vuqLdpYddXHDO73MyaRfH1\nM8xsspklmdkR9t2RYa8xs6lmVt/MHjGz46JVh8SXlLKbiITP3TcQGaYEM7sLyHf3P4VaVKSWZHff\nV8LTlxO5TuSbCrxein87zk1ZrgT+4e77i46DZmY/Aa4hcn3KFjN7DPgL3w4kJ1IibSlIzDOzS4Nx\n6Gea2ePBN+cUM9tsZn80sxlmNt7MjgoGCFxiwT01zOxKM3sjeH6Bmf2mnK97j5lNAfqZ2d3Bt/Iv\nzeyJ4Eri84iE2IH7PKSZ2YoiV6webWYTg+l7zOxvZjYBeC54j4eC955lZleW0PUL+fbq1gM1XwDc\nApzs7hsB3H0x0NzMmlTpf7zEJYWCxDQz60Zk6JL+HhlILIXI2PIQGQvmbXfvTWTgtLuAE4iMrvnb\nIi/TL1inN3CBRW5eU9brznD3fu7+CfCou/cFugfPDXb3V4lchXueu/d09wJK1ws4w90vBkYAa929\nH9CXyABwrQ/qdy2gpX93SI/DgYeIBMJavutzoH8ZNYho95HEvBOJ/OGcFuxCqc23QyjvdPcJwfRs\nImMA7TWz2UBOkdcY7+6bAMzsX0SGWEgp5XULgDeKrH+Cmd1KZNiBLCJDT7xVwX6McfddwfTJQGcz\nKxpC7YHlRdo35fuDn60BthG5x8BjBz23FjisgjVJAlIoSKwz4Fl3/9/vLIyMGln02/l+YHeR6aI/\n+weP9eJlvO5OD8aHMbM6RPbX93b3lWZ2D5FwKM5evt06P7jN9oP6dJ27v1PC6wDsLOE1hgAfmtna\nYGvlgFrBOiKl0u4jiXUTgXPNLAsKz1JqXcY6BzvZzBoEf+CHAR9V4HVrEwmZ9WZWl2/vBAaRb+11\ni8wvBfoE00XbHWw8cF0QQJhZx2DUzELuvg6odfCZV+6+BhhM5OY7JxZ5qgPfjg4sUiJtKUhMc/fZ\nFrnB+cRgOOU9RM68qchY+h8CLxO5EdNL7j4ToDyv6+4bzOwFIn9wl/Hdu8A9BzxtZjuJHLe4C3jK\nzL6h9GGd/wa0BmYGu67WEgmrg71D5DjB5INqWmxmPwL+bWbDiOw6yyFyXEGkVBolVRJacGZPN3e/\nKexaKsoit1m9zt1/Uka7c4Au7n539VQmsUy7j0RilLtPJXL8oKzfYwMeroaSJA5oS0FERAppS0FE\nRAopFEREpJBCQURECikURESkkEJBREQKKRRERKTQ/wMAlkrxT4hRqwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"num_points = 1000\n",
"Tarray = np.linspace(250, 300, num_points)\n",
"alpha = np.zeros_like(Tarray)\n",
"for n in range(num_points):\n",
" alpha[n] = albedo(Tarray[n])\n",
" \n",
"plt.plot(Tarray, alpha)\n",
"plt.xlabel('Temperature (K)')\n",
"plt.ylabel('albedo')\n",
"plt.grid()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0malbedo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mTarray\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m\u001b[0m in \u001b[0;36malbedo\u001b[0;34m(T, ao, ai, To, Ti)\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0malbedo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mao\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.289\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mai\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.7\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTo\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m293.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTi\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m260.\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;34m'''Given a single input temperature in Kelvin, return a single albedo'''\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0mT\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0mTi\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 9\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mai\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mTi\u001b[0m\u001b[0;34m<\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m<\u001b[0m\u001b[0mTo\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
]
}
],
"source": [
"albedo(Tarray)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYVOXZx/HvvRXYpS8sSFtEepGyoIJRiA0sEGPDbixY\nE43GRN/ktSRqjIktBmPs5Y0i0RhIxCAo2JUmgrSlCEiR3pa2lPv9Yw7ritvZ2bMz8/tc11ycc+Y5\nM/fD7O5vTnuOuTsiIiIASWEXICIiNYdCQURECikURESkkEJBREQKKRRERKSQQkFERAopFEREpJBC\nQURECikURESkUErYBVRUVlaW5+TkVGrd7du3k5GRUbUF1XDqc2JQnxPDofR5+vTp6929SVntYi4U\ncnJymDZtWqXWnTx5MgMHDqzagmo49TkxqM+J4VD6bGbLytNOu49ERKSQQkFERAopFEREpJBCQURE\nCikURESkUFRDwcwGm9kCM1tkZrcV8/zDZjYzeOSZ2eZo1iMiIqWL2impZpYMjAROAlYAU81srLvP\nPdDG3X9epP1PgV7RqkdERMoWzesU+gGL3H0JgJmNAoYBc0tofz5wZ7SKmbp0I/9cWMCMggXReosa\nKXXLXgaGXYSIxIxohkIL4Osi8yuAo4praGZtgLbAuyU8PwIYAZCdnc3kyZMrXMy4rwr49+ICWLyo\nwuvGKgcMp3HtdzksM3EOH+Xn51fqZySWqc+JoTr6HM1QsGKWeQlthwOvufu+4p509yeBJwFyc3O9\nMlf0DRwIpybYFZDr83dz7O8n8tHWBow8vXfY5VQbXemaGNTn6Ijm18cVQKsi8y2BVSW0HQ68EsVa\nElJWZjon5aTy5qzVzF21NexyRCQGRDMUpgLtzaytmaUR+cM/9uBGZtYRaAh8EsVaEtbgnFTq1Urh\noQmJdSxFRConaqHg7nuBG4DxwDxgtLvPMbPfmtnQIk3PB0a5e0m7luQQZKQaI447nInz1vL58k1h\nlyMiNVxUjz66+zh37+Du7dz93mDZHe4+tkibu9z9e9cwSNX5yYC2NMpI48G388IuRURquMQ5JSWB\nZaSncN3Adny4aD2fLN4QdjkiUoMpFBLERUe3IbteOg9NWID21IlISRQKCaJWajI3/LA9U5du4r28\ndWGXIyI1lEIhgZyX24qWDWvz4Nt52loQkWIpFBJIWkoSN57QntkrtzB+zpqwyxGRGkihkGDO7NWC\nw5tk8NCEBezbr60FEfkuhUKCSUlO4ucndiBvTT7/mVXSBeYikqgUCgnotO7N6dSsLg9PyGPPvv1h\nlyMiNYhCIQElJRm3nNyRpRt2MHra12WvICIJQ6GQoE7s3JTcNg15dOJCdhTsDbscEakhFAoJysy4\nbUgn1m7bzXMfLQ27HBGpIRQKCSw3pxEndm7KE5MXs2l7QdjliEgNoFBIcLee0on8gr08Pjlx7kgn\nIiVTKCS4js3qclbvlrzwyTJWbt4ZdjkiEjKFgvDzkzoA8MgEDa0tkugUCkKLBrW55Og2vD5jBXlr\ntoVdjoiESKEgAFw/6Agy0lJ44L+6badIIlMoCAANM9K4ZmA7Js5bw7SlG8MuR0RColCQQj8ZkEOT\nuun84b/zNbS2SIJSKEihOmkp3HhC5EY878xbG3Y5IhIChYJ8x3l9W9E2K4MHxs/X0NoiCUihIN+R\nmpzErad0JG9NPv/QYHkiCUehIN8zpFsz+rRpyIMT8ti+W4PliSQShYJ8j5nx69M6s27bbv72/pKw\nyxGRaqRQkGL1bt2Q03o058n3F/PNll1hlyMi1UShICW6bXAn9u+HB9/WBW0iiUKhICVq1agOlw3I\n4bUZK5izakvY5YhINVAoSKmuH3gE9Wunct+4ebqgTSQBKBSkVPXrpHLjCe35aNEGJi9YF3Y5IhJl\nCgUp04VHtSGncR3uHTePvfv2h12OiESRQkHKlJaSxG1DOrNobT6v6oI2kbimUJByOaVrNv1yGvHw\nhDy27doTdjkiEiUKBSmXAxe0rc8v4In3FoddjohEiUJByu3IVg0Y1vMwnv7gK1Zs2hF2OSISBQoF\nqZBfDu6EGfz+rflhlyIiUaBQkApp0aA21xzfjjdnreazJRvCLkdEqphCQSrs6uPacVj9Wtz977m6\n54JInFEoSIXVTkvm9lM7M3f1VkbrFFWRuKJQkEo5vUdz+uY05E/jF7Blp05RFYkXUQ0FMxtsZgvM\nbJGZ3VZCm3PNbK6ZzTGzl6NZj1QdM+POM7qycUcBj72zMOxyRKSKRC0UzCwZGAkMAboA55tZl4Pa\ntAduBwa4e1fgpmjVI1WvW4v6nJfbiuc/XsridflhlyMiVSCaWwr9gEXuvsTdC4BRwLCD2lwFjHT3\nTQDuvjaK9UgU3HJyR2qnJnPPf+aGXYqIVAGL1nDIZnY2MNjdrwzmLwaOcvcbirT5F5AHDACSgbvc\n/b/FvNYIYARAdnZ2n1GjRlWqpvz8fDIzMyu1bqyqjj6/9dUeXl1QwM/7pHNkk5Sovld56HNODOpz\nxQwaNGi6u+eW1S6av8FWzLKDEygFaA8MBFoCH5hZN3ff/J2V3J8EngTIzc31gQMHVqqgyZMnU9l1\nY1V19Ln/sfuZ8sj7jFkO1555HGkp4Z6/oM85MajP0RHN394VQKsi8y2BVcW0GePue9z9K2ABkZCQ\nGJKWksT/nt6ZJeu28+InS8MuR0QOQTRDYSrQ3szamlkaMBwYe1CbfwGDAMwsC+gALIliTRIlgzo2\nZWDHJjwycSFrt+4KuxwRqaSohYK77wVuAMYD84DR7j7HzH5rZkODZuOBDWY2F5gE3OruGjshBh04\nRbVg737uGzcv7HJEpJKielTQ3ccB4w5adkeRaQduDh4S49pmZXD18Yfz2LuLGN6vNUcf3jjskkSk\ngnRFs1Sp6wYeQYsGtblzzBz26NadIjFHoSBVqnZaMnee0YUFa7bxwsdLwy5HRCpIoSBV7qQu2QwK\nDjqv0UFnkZiiUJAqZ2bcNbQrBft00Fkk1igUJCraNM7gmuPbMWbmKj5ZrBPKRGKFQkGi5rqB7WjZ\nsDZ3jPlSB51FYoRCQaKmVmoyd53RlYVr83n+o6VhlyMi5aBQkKg6sUs2J3RqyiMT8/hmiw46i9R0\nCgWJujvP6Mre/c5v/zMn7FJEpAwKBYm61o3r8LMT2jNu9je8M29N2OWISCkUClItrvrB4XTIzuSO\nMXPYvntv2OWISAkUClIt0lKSuO/M7qzcvJNHJuaFXY6IlEChINUmN6cR5/drzbMfLeXLlVvCLkdE\niqFQkGp12+BONKyTxq/fmM2+/dG5FayIVJ5CQapV/Tqp3HFGF75YsYWXPlkadjkichCFglS7M3o0\n57gOTfjT23ms3rIz7HJEpAiFglQ7M+OeYd3Ys28/d4+dG3Y5IlKEQkFC0bpxHW48sT3/nfMNE+bq\n2gWRmkKhIKG56geH0zG7LneO+ZJ8XbsgUiMoFCQ0qclJ3Pfj7qzeuosH/js/7HJEBIWChKxPm4Zc\n1j+HFz9ZxpSvNoZdjkjCUyhI6G49pSOtGtXmV6/PYteefWGXI5LQFAoSujppKdz/4x58tX47D2sI\nDJFQKRSkRhhwRBbD+7biqfeXMGvF5rDLEUlYCgWpMf7ntM40qZvOL1+bRcFe3b5TJAwKBakx6tVK\n5d4fdWf+N9v46+TFYZcjkpAUClKjnNglm6FHHsZfJi1kwTfbwi5HJOEoFKTGufOMLtStlcovX5+l\nkVRFqplCQWqcxpnp3HlGF774ejNPf7Ak7HJEEopCQWqkoUcexsldsnnw7Tzy1mg3kkh1KXcomFma\nmXULHqnRLErEzLj3zO5k1krhltFfsGefzkYSqQ7lCgUzGwgsBEYCjwN5ZnZcFOsSoUnddO75UTdm\nr9zC45N0NpJIdSjvlsKDwMnufry7HwecAjwcvbJEIk7t3pxhPQ/jsXcX6r7OItWgvKGQ6u4LDsy4\nex6gXUhSLe4e2pVGGWncPHomu/dqbCSRaCpvKEwzs2fMbGDweAqYHs3CRA5oUCeNP5zdg7w1+Tw8\nYWHY5YjEtfKGwrXAHOBnwI3AXOCaaBUlcrBBHZtyfr9WPPn+YqYv0xDbItFSrlBw993u/pC7/9jd\nz3T3h919d7SLEynq16d14bAGtbll9BfsKNCd2kSiodRQMLPZZjarpEd1FSkCkJmewh/PPpKlG3bw\nh7d0pzaRaEgp4/nTg3+vD/59Kfj3QmBHVCoSKcUx7RrzkwE5PPfRUgZ1asrAjk3DLkkkrpS6peDu\ny9x9GTDA3X/p7rODx21ETkstlZkNNrMFZrbIzG4r5vnLzGydmc0MHldWviuSKH41uBMds+vyi3/M\nYkO+9mKKVKXyHmjOMLNjD8yYWX8go7QVzCyZyMVuQ4AuwPlm1qWYpq+6e8/g8XQ565EEVis1mUfP\n78nWXXv41euzcNegeSJVpbyhcAUw0syWmtlXRK5qvryMdfoBi9x9ibsXAKOAYZUvVeRbnZrV47bB\nnZg4by0vT1kedjkiccMq8i3LzOoF65R5aamZnQ0Mdvcrg/mLgaPc/YYibS4Dfg+sA/KAn7v718W8\n1ghgBEB2dnafUaNGlbvmovLz88nMzKzUurEqnvu8352Hpu0mb9M+7upfm8MyI99x4rnPJVGfE8Oh\n9HnQoEHT3T23zIbuXuYDyAaeAd4K5rsAV5SxzjnA00XmLwYeO6hNYyA9mL4GeLesWvr06eOVNWnS\npEqvG6vivc9rtuz0nneP91Mffd9379nn7vHf5+Koz4nhUPoMTPNy/L0v7+6j54HxwGHBfB5wUxnr\nrABaFZlvCaw6KJA2+LfXOzwF9ClnPSIANK1Xiz+c1YM5q7by4IQFZa8gIqUqbyhkuftoYD+Au+8F\nyhqEZirQ3szamlkaMBwYW7SBmTUvMjsUmFfOekQKndy1GRcc1Zon31/Cx4vXh12OSEwrbyhsN7PG\ngAOY2dFAqccVguC4gcgWxjxgtLvPMbPfmtnQoNnPzGyOmX1BZAiNyyrRBxF+c1pn2mZlcPOrX5Bf\noLORRCqrrIvXDriZyLf8dmb2EdAEOLusldx9HDDuoGV3FJm+Hbi93NWKlKBOWgp/Ht6LMx//iGe+\nTOK0kxwzC7sskZhT3rGPZgDHA/2Bq4Gu7q5hLqRG6daiPrcN6czna/fx/MdLwy5HJCaV985rtYjs\n3vkdcDdwfbBMpEa5fEAOvZomc9+4ecxasTnsckRiTnmPKbwIdAUeA/5C5JTUl0pdQyQEZsYV3dJp\nkpnODS9/ztZde8IuSSSmlDcUOrr7Fe4+KXiMADpEszCRyspMMx67oBcrN+/k9n/O1jAYIhVQ3lD4\nPDjjCAAzOwr4KDoliRy6Pm0a8YuTO/LmrNUaBkOkAko9+8jMZhM5DTUVuMTMlgfzbYjcfU2kxrr6\nuMP5ZMkG7v73XHq3bkjn5vXCLkmkxitrS+F04AxgMNCWyBlIA4Pp06JamcghSkoyHjr3SBrUTuX6\nl2ewfbfu1iZSlrJCYVsZD5EaLSsznUeH92Lp+u06viBSDmVdvDadyO6iA1cBHfiNsmD68CjVJVJl\njmnXmFtO7sgfxy+gd+sGXDagbdglidRYpYaCuxf+9phZI6A9oOsTJOZce3w7Pl++iXvenEf3lvXp\n06ZR2CWJ1EjlvXjtSuA94L/AXcG/d5S2jkhNkpRkPHhuT1o0rM11f5/Bum26jadIccp7SuqNQF9g\nmbsPAnoBGo5SYkr92qn89cI+bNm5h5++MoO9+/aHXZJIjVPeUNjl7rsAzCzd3ecDHaNXlkh0dDms\nHvf+qDufLtnIH9/W/RdEDlbeUVJXmFkD4F/ABDPbxEE3zBGJFWf1acmM5Zv423tL6NWqIYO7NQu7\nJJEao1yh4O5nBpN3mdkkoD6R4woiMemOM7rw5aqt3PqPL+iQncnhTRLrXr8iJSnv7qNC7v6eu491\n94JoFCRSHdJTknn8wt6kpiQx4qXpbNPAeSJAJUJBJF60aFCbkRf0Zun67dw0aib79+vCNhGFgiS0\nY9o15s4zuvDO/LU8OEEHnkXKe6BZJG5ddHQb5q7eyshJi+nUrB5nHHlY2CWJhEZbCpLwzIy7h3Yj\nt01Dbn3tC75cuSXskkRCo1AQAdJSkvjrRX1oWCeNES9OY32+rniWxKRQEAk0qZvOkxfnsmF7Adf9\n3wwK9uqKZ0k8CgWRIrq3rM8DZ/dgytKN3DHmSw21LQlHB5pFDjKsZwvy1mxj5KTF5GRlcM3x7cIu\nSaTaKBREinHLSR1ZtmEH9781nzaN6jCke/OwSxKpFtp9JFKMpCTjT+ccSe/WDbjp1Zl8vnxT2CWJ\nVAuFgkgJaqUm89QluTStl85VL07j6407wi5JJOoUCiKlaJyZznOX9WX33v1c/vxUtuzUGEkS3xQK\nImU4omld/nZRH75av53r/z6DPbo5j8QxhYJIOfQ/Iov7zuzOh4vW85s3dKqqxC+dfSRSTuf2bcXy\njTv4y6RFZNdL5+aTdfNBiT8KBZEKuOXkDqzdtos/v7uIJvVqcfHRbcIuSaRKKRREKsDMuO/M7mzI\nL+COMV+SlZGmaxgkruiYgkgFpSQn8ZcLetOrVQNuHDWTT5dsCLskkSqjUBCphNppyTxzaV9aN67D\nVS9MY97qrWGXJFIlFAoildQwI40XLu9HRnoKlz47RRe3SVxQKIgcghYNavPC5f3YtWcfFz3zGWu3\n7gq7JJFDolAQOUQdm9Xl+cv7sW7bbi58+jM2bi8IuySRSlMoiFSB3q0b8vSluSzfuINLnv2Mrbs0\nHIbEpqiGgpkNNrMFZrbIzG4rpd3ZZuZmlhvNekSiqX+7LP56UW/mr97G5c9NZUfB3rBLEqmwqIWC\nmSUDI4EhQBfgfDPrUky7usDPgM+iVYtIdflhp2weHd6LGcs3cfVL09m1Z1/YJYlUSDS3FPoBi9x9\nibsXAKOAYcW0+x3wAKAjdBIXTuvRnD+c1YMPFq7np698rgH0JKZEMxRaAF8XmV8RLCtkZr2AVu7+\nnyjWIVLtzsltxd1DuzJh7hp+pmCQGGLRGu3RzM4BTnH3K4P5i4F+7v7TYD4JeBe4zN2Xmtlk4Bfu\nPq2Y1xoBjADIzs7uM2rUqErVlJ+fT2ZmZqXWjVXqc7jGL93DK/ML6JOdzLVHppOSZFF5n5rU5+qi\nPlfMoEGDprt72cdt3T0qD+AYYHyR+duB24vM1wfWA0uDxy5gFZBb2uv26dPHK2vSpEmVXjdWqc/h\ne/qDJd7mV//xES9O9d179kXlPWpan6uD+lwxwDQvx9/uaO4+mgq0N7O2ZpYGDAfGFgmjLe6e5e45\n7p4DfAoM9WK2FERi2RXHtuWO07swfs4abnh5BgV7tStJaq6ohYK77wVuAMYD84DR7j7HzH5rZkOj\n9b4iNdHlx7blzjO68PbcNVyvYJAaLKpDZ7v7OGDcQcvuKKHtwGjWIhK2nwxoS5IZd46dw3V/n8HI\nC3uRnpIcdlki36ErmkWq0aX9c/jdsK5MnLeGK1+YpgvcpMZRKIhUs4uPyeGPZ/fgo0XrufiZKWzZ\nqSExpOZQKIiE4JzcVoy8oDezVmzm/Cc/ZX3+7rBLEgEUCiKhGdK9OU9f2pcl6/M592+fsGrzzrBL\nElEoiITp+A5NeOmKo1i3dTfnPPEJS9blh12SJDiFgkjI+uY04pURR7Nzzz7OfuITPl++KeySJIEp\nFERqgG4t6vP6tf3JTE/h/Kc+5Z15a8IuSRKUQkGkhmiblcHr1/anQ3ZdrnpxGi9/tjzskiQBKRRE\napAmddN55aqjOb5DE/7njdk89PaCA2OFiVQLhYJIDZORnsJTl+RyXm4r/vzuIm59bZaGxZBqE9Vh\nLkSkclKSk7j/rO40b1CLRyYuZPnGHTxxUR8aZaSFXZrEOW0piNRQZsZNJ3bg0eE9mfn1ZoaN/JC8\nNdvCLkvinEJBpIYb1rMFo68+hl179vPjxz/m3fk6M0miR6EgEgN6tmrA2BsGkJNVhytemMZT7y/R\nAWiJCoWCSIxoXr82o68+hiHdmnHvuHncMvoLdhbsC7ssiTMKBZEYUicthb+c35ufn9iBN2au5MzH\nP2LZhu1hlyVxRKEgEmOSkowbT2zPs5f1ZfWWXZz+2IfMXKv7MkjVUCiIxKhBHZvyn58eS+tGdXhk\nxm4efHsB+/brOIMcGoWCSAxr1agOr1/bnx+0SOGxdxdx2XNTWLdN92aQylMoiMS4WqnJXN4tjd//\nuDtTvtrIkEc/4P28dWGXJTFKoSASB8yM8/u1ZuwNx9IoI5VLnp3C78fN0/AYUmEKBZE40rFZXcZc\nfywXHtWav72/hLOf+Jil63V2kpSfQkEkztROS+beM7vzxEW9Wbp+O6f9+QNGTVmui92kXBQKInFq\ncLfmvHXTcfRo2YDb/jmbnzw/lW+27Aq7LKnhFAoicaxFg9r8/cqjuHtoVz5bspGTH36Pf85Yoa0G\nKZFCQSTOJSUZl/bP4a0bf0CH7LrcPPoLrn5pOmu3aqtBvk+hIJIgcrIyePXqY/j1qZ2ZnLeOEx58\nj5c+XcZ+XfAmRSgURBJIcpJx1XGHM/6m4+jRqj7/+68vOeuJj5n/zdawS5MaQqEgkoDaZmXwf1cc\nxcPnHcmyDTs4/c8fcv9b8zXqqigURBKVmXFmr5a8c/Px/Lh3C554bzEnPDiZMTNX6kB0AlMoiCS4\nhhlpPHD2kfzjmmNolJnGjaNmcvYTn/DF15vDLk1CoFAQEQD65jRizPXH8sBZPVi2YTvDRn7ELaO/\nYI3OUkooKWEXICI1R3KScW7fVgzp3oy/TFrEcx8u5c3Zq7i0fw7XHNeOhhlpYZcoUaYtBRH5nrq1\nUrl9SGcm3nw8Q7o158n3l3DcA5N47J2F5O/WDX3imUJBRErUunEdHj6vJ/+98TiOadeYByfkcfwD\nk3j6gyXsKFA4xCOFgoiUqWOzujx5SS5vXNefTs3rcs+b8xhw/7s8OnEhm3cUhF2eVCGFgoiUW6/W\nDfn7lUfz+rX96dOmIQ9PzGPA/e9y75tzdUA6TuhAs4hUWJ82DXn60r7M/2Yrf528mGc+/IrnP17K\nqd2bc2n/HHq1aoCZhV2mVIJCQUQqrVOzejw6vBe3nNSRZz/6itemr2DMzFX0aFmfS47J4fQezamV\nmhx2mVIBUd19ZGaDzWyBmS0ys9uKef4aM5ttZjPN7EMz6xLNekQkOlo3rsNdQ7vy6f+cwO+GdWVH\nwT5+8Y8vOOb373DX2Dl8uXKLrpKOEVHbUjCzZGAkcBKwAphqZmPdfW6RZi+7+xNB+6HAQ8DgaNUk\nItGVmZ7CxcfkcNHRbfh48QZembKcl6cs5/mPl9KpWV3O7tOSYT1b0KRuetilSgmiufuoH7DI3ZcA\nmNkoYBhQGAruXnRoxgxAXyVE4oCZMeCILAYckcWWHXsYO2sVr01fwT1vzuO+cfPo17YRp3ZvzuCu\nzWhar1bY5UoR0QyFFsDXReZXAEcd3MjMrgduBtKAH0axHhEJQf06qVx8dBsuProNC9ds49+zVjNu\n9mruGDOHO8fOIbdNQ07u0ozjOzahfdNMHaAOmUVrP5+ZnQOc4u5XBvMXA/3c/acltL8gaH9pMc+N\nAEYAZGdn9xk1alSlasrPzyczM7NS68Yq9TkxxGKfV+bvZ9o3e5n6zV5W5Ef+DjWqZXTLSqZ7VjKd\nGyWTmVZyQMRinw/VofR50KBB0909t6x20QyFY4C73P2UYP52AHf/fQntk4BN7l6/tNfNzc31adOm\nVaqmyZMnM3DgwEqtG6vU58QQ631euXkn7+et470F6/ho0Xq2BUNptG+aSW5OQ3LbNCI3pyGtG9Up\n3JKI9T5XxqH02czKFQrR3H00FWhvZm2BlcBw4IKiDcysvbsvDGZPAxYiIgmnRYPanN+vNef3a82e\nffuZ+fVmpny1kWlLN/KfWat5ZUpkT3S9Wil0OawenZvXI2nLHrJWbqFtVgYZ6Tq7vqpE7X/S3fea\n2Q3AeCAZeNbd55jZb4Fp7j4WuMHMTgT2AJuA7+06EpHEkpqcRN+cRvTNaQTA/v1O3tptTF+2iTmr\ntjJ31VZGTfmanXv28cyXHwKQlZlG60Z1yGmcQatGdWhSN52szHSa1E0jKzOdRhlp1ElLITlJxyvK\nEtV4dfdxwLiDlt1RZPrGaL6/iMS+pCSjU7N6dGpWr3DZvv3O6HGTqNe6M0s3bGf5hh0s27idT5ds\n4I2ZKylpr3haShJ10pKpk5pMWkoSZkZhTBiF099ZXoOc2HwvA6P8HtrmEpGYk5xkNM9MYmCP5t97\nbs++/WzIL2B9/m7W5e9m/bbdbNpRwM6C/ezYs5edBfvYUbCPgr37C8+Bd/dvz4d38Bp6dnxG6qao\nv4dCQUTiSmpyEs3q16JZ/fi7/mHy5MlRfw+NkioiIoUUCiIiUkihICIihRQKIiJSSKEgIiKFFAoi\nIlJIoSAiIoUUCiIiUihqo6RGi5mtA5ZVcvUsYH0VlhML1OfEoD4nhkPpcxt3b1JWo5gLhUNhZtPK\nM3RsPFGfE4P6nBiqo8/afSQiIoUUCiIiUijRQuHJsAsIgfqcGNTnxBD1PifUMQURESldom0piIhI\nKeIqFMyslZlNMrN5ZjbHzG4Mlt9lZivNbGbwOLXIOreb2SIzW2Bmp4RXfcWV1N/guZ8GfZpjZg8U\nWR6z/YVSP+NXi3y+S81sZpF14rXPPc3s06DP08ysX7DczOzPQZ9nmVnvcHtQcaX0+Ugz+8TMZpvZ\nv82sXpF1Yv1zrmVmU8zsi6DPdwfL25rZZ2a2MPg5TwuWpwfzi4Lnc6qkEHePmwfQHOgdTNcF8oAu\nwF3AL4pp3wX4AkgH2gKLgeSw+1EF/R0ETATSg+eaxkN/S+vzQW0eBO6I9z4DbwNDguWnApOLTL9F\n5O6SRwOfhd2HKuzzVOD4YPnlwO/i6HM2IDOYTgU+Cz6/0cDwYPkTwLXB9HXAE8H0cODVqqgjrrYU\n3H21u88IprcB84AWpawyDBjl7rvd/StgEdAv+pVWjVL6ey1wv7vvDp5bG6wS0/2Fsj9jMzPgXOCV\nYFE899kYvA6kAAAGiElEQVSBA9+U6wOrgulhwIse8SnQwMy+f9/KGqyUPncE3g+aTQDOCqbj4XN2\nd88PZlODhwM/BF4Llr8A/CiYHhbMEzx/QvDzf0jiKhSKCjalehFJW4Abgk3pZ82sYbCsBfB1kdVW\nUHqI1FgH9bcD8INgk/I9M+sbNIub/kKxnzHAD4A17r4wmI/nPt8E/NHMvgb+BNweNIvnPn8JDA2e\nOgdoFUzHRZ/NLDnY9bmWSOgtBja7+96gSdF+FfY5eH4L0PhQa4jLUDCzTOB14CZ33wr8FWgH9ARW\nE9m9AJHNtYPF3OlYxfQ3BWhIZNPzVmB08A0iLvoLxfb5gPP5disB4rvP1wI/d/dWwM+BZw40LWb1\neOnz5cD1ZjadyG6lggNNi1k95vrs7vvcvSfQksiWTufimgX/RqXPcRcKZpZK5Ifo7+7+TwB3XxP8\nZ+8HnuLbzcoVfPtNAyIfxCpiSHH9JdKvfwabo1OA/UTGTIn5/kKJfcbMUoAfA68WaR7Pfb4UODD9\nD+Lo5xpK/F2e7+4nu3sfIuG/OGgeF30+wN03A5OJfLFrEPxsw3f7Vdjn4Pn6wMZDfe+4CoXg2/Az\nwDx3f6jI8qL7U88ksgkKMBYYHhzFbwu0B6ZUV72HqqT+Av8ish8SM+sApBEZRCum+wul9hngRGC+\nu68osiye+7wKOD6Y/iFwYJfZWOCS4Cyko4Et7r662gquAqX8LjcN/k0CfkPkwCvEx+fcxMwaBNO1\nifw8zwMmAWcHzS4FxgTTY4N5guff9eCo8yEJ+4h7VT6AY4lsPs0CZgaPU4GXgNnB8rFA8yLr/JrI\nt40FBGdyxMqjlP6mAf9HJPxmAD+Mh/6W1ufgueeBa4pZJy77HCyfTuSsm8+APkF7A0YGfZ4N5Ibd\nhyrs841EzkTKA+4nuAA3Tj7nHsDnQZ+/5Nsz6A4nEnCLiGwRHjirsFYwvyh4/vCqqENXNIuISKG4\n2n0kIiKHRqEgIiKFFAoiIlJIoSAiIoUUCiIiUkihIDHBzBoXGQX1G/vuqLdpYddXHDO73MyaRfH1\nM8xsspklmdkR9t2RYa8xs6lmVt/MHjGz46JVh8SXlLKbiITP3TcQGaYEM7sLyHf3P4VaVKSWZHff\nV8LTlxO5TuSbCrxein87zk1ZrgT+4e77i46DZmY/Aa4hcn3KFjN7DPgL3w4kJ1IibSlIzDOzS4Nx\n6Gea2ePBN+cUM9tsZn80sxlmNt7MjgoGCFxiwT01zOxKM3sjeH6Bmf2mnK97j5lNAfqZ2d3Bt/Iv\nzeyJ4Eri84iE2IH7PKSZ2YoiV6webWYTg+l7zOxvZjYBeC54j4eC955lZleW0PUL+fbq1gM1XwDc\nApzs7hsB3H0x0NzMmlTpf7zEJYWCxDQz60Zk6JL+HhlILIXI2PIQGQvmbXfvTWTgtLuAE4iMrvnb\nIi/TL1inN3CBRW5eU9brznD3fu7+CfCou/cFugfPDXb3V4lchXueu/d09wJK1ws4w90vBkYAa929\nH9CXyABwrQ/qdy2gpX93SI/DgYeIBMJavutzoH8ZNYho95HEvBOJ/OGcFuxCqc23QyjvdPcJwfRs\nImMA7TWz2UBOkdcY7+6bAMzsX0SGWEgp5XULgDeKrH+Cmd1KZNiBLCJDT7xVwX6McfddwfTJQGcz\nKxpC7YHlRdo35fuDn60BthG5x8BjBz23FjisgjVJAlIoSKwz4Fl3/9/vLIyMGln02/l+YHeR6aI/\n+weP9eJlvO5OD8aHMbM6RPbX93b3lWZ2D5FwKM5evt06P7jN9oP6dJ27v1PC6wDsLOE1hgAfmtna\nYGvlgFrBOiKl0u4jiXUTgXPNLAsKz1JqXcY6BzvZzBoEf+CHAR9V4HVrEwmZ9WZWl2/vBAaRb+11\ni8wvBfoE00XbHWw8cF0QQJhZx2DUzELuvg6odfCZV+6+BhhM5OY7JxZ5qgPfjg4sUiJtKUhMc/fZ\nFrnB+cRgOOU9RM68qchY+h8CLxO5EdNL7j4ToDyv6+4bzOwFIn9wl/Hdu8A9BzxtZjuJHLe4C3jK\nzL6h9GGd/wa0BmYGu67WEgmrg71D5DjB5INqWmxmPwL+bWbDiOw6yyFyXEGkVBolVRJacGZPN3e/\nKexaKsoit1m9zt1/Uka7c4Au7n539VQmsUy7j0RilLtPJXL8oKzfYwMeroaSJA5oS0FERAppS0FE\nRAopFEREpJBCQURECikURESkkEJBREQKKRRERKTQ/wMAlkrxT4hRqwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# A vectorized version using numpy.where (look it up!)\n",
"\n",
"def albedo_vectorized(T, ao=0.289, ai=0.7, To=293., Ti=260.):\n",
" alpha1 = np.where(T"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(Years, Tsteps_C1, color='red', label='C1 (100 meters)')\n",
"plt.plot(Years, Tsteps_C2, color='blue', label='C2 (200 meters)')\n",
"plt.plot(Years, Tsteps_varalpha, color='magenta', label='C1, variable alpha')\n",
"plt.xlabel('Years')\n",
"plt.ylabel('Global mean temperature (K)');\n",
"plt.legend()\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bonus problem\n",
"\n",
"*Open-ended investigation for extra credit, not required*\n",
"\n",
"Something very different occurs in this model if you introduce a strong negative radiative forcing, either by substantially reducing greenhouse gases (which we would represent as an increase in the transmissivity $\\tau$), or by decreasing the incoming solar radiation $Q$.\n",
"\n",
"Investigate, using your numerical model code, and report your results along with your thoughts."
]
},
{
"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.6.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}