<jupyter_start><jupyter_text>Introducing the Column Model  About this exercise: There are two goals for this (hopefully straightforward) exercise:1. To introduce you to working with files and modules in Python2. To introduce the column radiation modelWe will be working with this same model code for the next few weeks.Your job is to go through the whole exercise. There are 3 questions clearly marked below for you to answer. Hand in ONLY your answers to those questions.  Please DO NOT just copy and paste the examples from here. Type them in yourself! This is MUCH better way to learn what you are doing.  Answers are due on Tuesday February 25, 2014  About Python modules: Every time we use the import statement in Python, the interpreter loads instructions from a file, usually with the extension .py at the end of the file name.For example if we typed 'import foo' the interpreter would try to open and read a file called 'foo.py'This file would contain the same types of Python code that you have been entering manually on the command line, e.g. definitions of variables and functions.The real beauty of the import statement is that it makes the Python language very modular. With one 'import' we can have access to all kinds of useful code.To then access variables and functions in a module, we usually type the name of the module, followed by period (.), followed by the variable or function name. You've already seen a bit of this with the netCDF4 module you used in the last homework. You'll practice these skills again here.<jupyter_code>%matplotlib inline
import numpy as np
import climlab
from climlab import constants as const<jupyter_output><empty_output><jupyter_text>Creating a new instance of the Column Model To begin working with the model, we just need to tell Python to create a new object called a `mycolumn` that implements a grey radiation model.Try this:<jupyter_code>mycolumn = climlab.GreyRadiationModel( num_lev=2 )<jupyter_output><empty_output><jupyter_text>You have just created a new column object. What is that? Let's take a look:<jupyter_code>print(mycolumn)<jupyter_output>climlab Process of type <class 'climlab.model.column.GreyRadiationModel'>. 
State variables and domain shapes: 
  Ts: (1,) 
  Tatm: (2,) 
The subprocess tree: 
Untitled: <class 'climlab.model.column.GreyRadiationModel'>
   LW: <class 'climlab.radiation.greygas.GreyGas'>
   SW: <class 'climlab.radiation.greygas.GreyGasSW'>
   insolation: <class 'climlab.radiation.insolation.FixedInsolation'><jupyter_text>Try just typing mycolumn. and then hitting the tab key. You should see a list pop up. Everything on that list is the name of some object that is contained within mycolumn."Object" is a very general concept. Everything in Python that has a name is an object. This includes variables (numbers, character strings, etc.), but also functions.Our object called mycolumn contains a bunch of data, and it also contains some functions that operate on those data.Try this:<jupyter_code>print(mycolumn.state)<jupyter_output>{'Ts': Field([288.]), 'Tatm': Field([200., 278.])}<jupyter_text>We will use the column object to model temperature in a column. What you just did above was to look at the temperatures currently stored in mycolumn, for the surface (a single number), and the atmosphere.When you created mycolumn a few lines ago, you specifically asked for a column with 2 levels in the atmosphere.We can see where those levels are in the atmosphere as follows:<jupyter_code>#  Three ways to get this:
print(mycolumn.lev)
print(mycolumn.Tatm.domain.axes['lev'].points)
print(mycolumn.domains['Tatm'].axes['lev'].points)<jupyter_output>[250. 750.]
[250. 750.]
[250. 750.]<jupyter_text>which is actually the pressure (in mb or hPa) at which we are specifying the temperature -- in this case essentially the lower and upper troposphere.  So what does this code actually do? It calculates the warming and cooling of the air and the surface based on the grey radition approximation that we have looked at in class.I encourage you to look through the code. Try typing this:<jupyter_code>edit ColumnModel<jupyter_output><empty_output><jupyter_text>You should find that the code is now displayed in your editor window. Have a look, but obviously don't fret about not understanding it at this point.Notice, though, that there are lots of comments sprinkled throughout the code. Comments in Python begin with  and are just words to help us understand the code. Comments are never executed by the Python interpreter.Convince yourself of this by typing something like this:<jupyter_code># this is not valid Python code<jupyter_output><empty_output><jupyter_text>You might want to try typing the same thing without the  in front.<jupyter_code>this is not valid Python code<jupyter_output><empty_output><jupyter_text>OK but what does this code actually do? It calculates the warming and cooling of the air and the surface based on the grey radition approximation that we have looked at in class.Try this:<jupyter_code>mycolumn.compute_diagnostics()
print(mycolumn.LW_absorbed_sfc)
print(mycolumn.LW_absorbed_atm)<jupyter_output>[-205.79241571]
[  87.88678511 -116.45384996]<jupyter_text>What you just did was to call a function Longwave_Heating() that calculates how much longwave radiation is emitted up and down from each layer, and how much is absorbed by each layer (and the surface).You then printed out some quantities that were calculated and stored by the function, which are actually the heating rates at each level in W / m$^2$Now try this:<jupyter_code>print(mycolumn.SW_absorbed_sfc)
print(mycolumn.SW_absorbed_atm)<jupyter_output>[239.2513]
[0. 0.]<jupyter_text>Hopefully this makes sense. Our code also calculates the energy absorbed due to shortwave (solar) radiation. The atmosphere is transparent to solar radiation in this model, so the absorption all occurs at the surface.Look in the code file to find where the function Shortwave_Heating() is defined. It's really just calculating$(1-\alpha) Q$ based on some default parameter values for $\alpha$ and $Q$. We'll think about changing the defaults later.  Energy (im)balance in the column To get the total energy sources and sinks at each point, we just need to add up the shortwave and longwave terms:<jupyter_code>print(mycolumn.LW_absorbed_atm + mycolumn.SW_absorbed_atm)
print(mycolumn.LW_absorbed_sfc + mycolumn.SW_absorbed_sfc)<jupyter_output>[  87.88678511 -116.45384996]
[33.45888429]<jupyter_text>Evidently this column is NOT in energy balance! The surface is gaining energy at a rate 33 W / m$^2$, the lower atmosphere is losing energy at 116 W / m$^2$, and the upper atmosphere is gaining nearly 90 W / m$^2$.  OK so what? The code will not just calculate the energy imbalance, but also change the temperatures in response to the imbalance. It does this by time-stepping, just like we did with the zero-dimensional model in the first homework.<jupyter_code>mycolumn.step_forward()
print(mycolumn.state)<jupyter_output>{'Ts': Field([288.69137531]), 'Tatm': Field([201.48238045, 276.0357808 ])}<jupyter_text>Here you just called a function `step_forward()` that computes the energy imbalance as you just did above, and then uses that imbalance to adjust the temperatures.  HOMEWORK QUESTION 1 Have the temperatures gone up or down at the surface and at each level from where they started? Why?(This is an easy question, not a trick question)  Timestepping to equilibrium Just like we did with the zero-dimensional model in the first homework, we will use loops to time-step the model towards equilibrium. Try this:<jupyter_code>for n in range(10):
    mycolumn.step_forward()
    print(mycolumn.Ts)<jupyter_output>[289.22566617]
[289.63070626]
[289.92977194]
[290.14229416]
[290.28446635]
[290.36976133]
[290.40936944]
[290.41256879]
[290.38703728]
[290.33911515]<jupyter_text>This little loop just repeated the call to Step_Forward 10 times, and printed out the surface temperature after each time step.Notice that the temperature is changing each time. That means we are not at equilibrium. Try it again!<jupyter_code>for n in range(10):
    mycolumn.step_forward()
    print(mycolumn.Ts)<jupyter_output>[290.27402541]
[290.19605881]
[290.10872876]
[290.01490115]
[289.91690295]
[289.81661304]
[289.71553825]
[289.61487688]
[289.51557179]
[289.41835483]<jupyter_text>Still changing, but not by as much.Here's a trick:<jupyter_code>mycolumn.integrate_years(0.3)
print(mycolumn.Ts)<jupyter_output>Integrating for 109 steps, 109.57266 days, or 0.3 years.
Total elapsed time is 0.35592820325800245 years.
[287.80463186]<jupyter_text>What you just did was to loop through the time-stepping procedure 100 times!Look at the code and find the function Step_Forward(). Do you see how the code creates a loop?In this case the function Step_Forward() takes an optional input argument which is the number of iterations through the loop. This number defaults to 1 if we don't specify it, which is what happened above!Has the model reached equilibrium yet? We can always keep on time-stepping and see if anything changes:<jupyter_code>mycolumn.integrate_years(0.3)
print(mycolumn.Ts)<jupyter_output>Integrating for 109 steps, 109.57266 days, or 0.3 years.
Total elapsed time is 0.6543603121435584 years.
[287.84175315]<jupyter_text>The changes are now minimal, and it is close to equilibrium.Let's look at the whole column temperature:<jupyter_code>for key, item in mycolumn.state.items():
    print(key, item)
mycolumn.diagnostics<jupyter_output><empty_output><jupyter_text>HOMEWORK QUESTION 2 Compare the temperatures you found here (after time-stepping to equilibrium) with the radiative equilibrium temperatures we derived in class for this same model. Do they agree?  Greenhouse warming in the 2-layer model Now that our column is in equilibrium, let's look at the Outgoing Longwave Radiation. The model keeps track of this for us:<jupyter_code>print(mycolumn.OLR)<jupyter_output>[239.22616887]<jupyter_text>This should hopefully be almost exactly equal to the shortwave absorption:<jupyter_code>print(mycolumn.ASR)<jupyter_output>[239.2513]<jupyter_text>Now you are going to do a "global warming" experiment, like we started in class.The following will increase the emissivity / absorptivity of each layer by 10%, which is analagous to an increase in greenhouse gases in the atmosphere:<jupyter_code>#  Make an exact clone of the model
column2 = climlab.process_like(mycolumn)
absorptivity = column2.subprocess['LW'].absorptivity
print(absorptivity)
column2.subprocess['LW'].absorptivity *= 1.1 
print(column2.subprocess['LW'].absorptivity)<jupyter_output>[0.47737425 0.47737425]
[0.52511167 0.52511167]<jupyter_text>Let's now re-calculate the longwave radiation with this new value of eps:<jupyter_code>column2.compute_diagnostics()<jupyter_output><empty_output><jupyter_text>HOMEWORK QUESTION 3 Find the new value of OLR after this change. Is it larger or smaller than it was before we added greenhouse gases? What do you think should happen to the surface temperature as a result? Why?<jupyter_code>print(column2.OLR)
column2.diagnostics<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06859400251391633, 'fraction_numerical': 0.04875202011133058, 'mean_word_length': 5.5872264931992905, 'pattern_counts': {'":': 0, '<': 80, '<?xml version=': 0, '>': 80, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '47301', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3615, 'n_tokens_neox': 3164, 'n_words': 1411, 'path': 'courseware/Introducing the Column Model.ipynb', 'repo_name': 'nfeldl/climlab', 'size': '11138'}