<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Using climlab to solve the 1D diffusive Energy Balance Model
An interactive tutorial for ENV / ATM 415
Step1:  What is a climlab Process?
climlab is a flexible Python engine for process-oriented climate modeling.
Every model in climlab consists of some state variables, and processes that operate on those state variables.
The state variables are the climate variables that we will step forward in time. For our EBM, it is just the surface temperature
Step2:  We have an array of temperatures in degrees Celsius.  Let's see how big this array is
Step3:  Every state variable exists on a spatial grid. In this case, the grid is an array of latitude points
Step4:  We can, for example, plot the current temperature versus latitude
Step5:  It is based on a very general concept of a model as a collection of individual, interacting processes. 
Every model in climlab is a collection of individual, interacting physical processes. 
In this case, the EBM_annual object has a list of 4 sub-processes. Each one is basically one of the terms in our energy budget equation above
Step6:  So what does it do?
Well for one thing, we can step the model forward in time. This is basically just like we did manually with our zero-dimensional model
Step7:  Looks like the temperature got a bit colder near the equator and warmer near the poles
Step8:  This value is in seconds. It is actually 1/90th of a year (so, to step forward one year, we need 90 individual steps). This is a default value -- we could change it if we wanted to.
Actually we can see a bunch of important parameters for this model all together in a dictionary called param
Step9:  Accessing the model diagnostics
Each model can have lots of different diagnostic quantities. This means anything that can be calculated from the current model state -- in our case, from the temperature.
For the EBM, this includes (among other things) the OLR and the ASR.
Step10:  This plot shows that $ASR > OLR$ (system is gaining extra energy) across the tropics, and $ASR < OLR$ (system is losing energy) near the poles.
It's interesting to ask how close this model is to global energy balance.
For that, we have to take an area-weighted global average.
Step11:  Running the model out to equilibrium
We write a loop to integrate the model through many timesteps, and then check again for energy balance
Step12:  Since there is still a significant imbalance, we are not yet at equilibrium. We should step forward again.
This time, let's use a convenient built-in function instead of writing our own loop
Step13:  We are now quite close to equilibrium. Let's make some plots
Step14:  Let's make some nice plots of all the terms in the energy budget.
We'll define a reusable function to make the plots.
Step15:  What if we start from a very different initial temperature?
Step16:  Why is the heat transport zero everywhere?
Now let's look at an animation showing how all the terms in the energy budget evolve from this initial condition toward equilibrium.
The code below generates a series of snapshots, that we can stitch together into an animation with a graphics program.
<ASSISTANT_TASK:>
Python Code:
#  We start with the usual import statements
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
#  create a new model with all default parameters (except the grid size)
mymodel = climlab.EBM_annual(num_lat = 30)
#  What did we just do?
print mymodel
Explanation: Using climlab to solve the 1D diffusive Energy Balance Model
An interactive tutorial for ENV / ATM 415: Climate Laboratory
Spring 2016
The model
The equation for the model is
$$ C \frac{\partial T}{\partial t} = ASR - OLR + \frac{D}{\cos\phi} \frac{\partial}{\partial \phi} \left( \cos\phi \frac{\partial T}{\partial \phi} \right)$$
Where, following our standard notation:
- $T$ is the surface temperature (in this case, the zonal average temperature)
- $C$ is the heat capacity in J m$^{-2}$ K$^{-1}$
- $ASR$ is the absorbed solar radiation
- $OLR$ is the outgoing longwave radiation
- $D$ is the thermal diffusivity in units of W m$^{-2}$ K$^{-1}$
- $\phi$ is latitude
Longwave parameterization
We will use a very simple parameterization linking surface temperature to OLR:
$$ OLR = A + B~T $$
Shortwave parameterization
$$ ASR = (1-\alpha) Q $$
and we will use the annual average insolation (ignore seasons).
The albedo $\alpha$ (at least for now) will be a simple prescribed function of latitude (larger at the poles than at the equator) that mimics the observed annual mean albedo.
The EBM is already coded up in climlab.
End of explanation
mymodel.Ts
Explanation: What is a climlab Process?
climlab is a flexible Python engine for process-oriented climate modeling.
Every model in climlab consists of some state variables, and processes that operate on those state variables.
The state variables are the climate variables that we will step forward in time. For our EBM, it is just the surface temperature:
End of explanation
mymodel.Ts.shape
Explanation: We have an array of temperatures in degrees Celsius.  Let's see how big this array is:
End of explanation
mymodel.lat
Explanation: Every state variable exists on a spatial grid. In this case, the grid is an array of latitude points:
End of explanation
plt.plot(mymodel.lat, mymodel.Ts)
plt.xlabel('Latitude')
plt.ylabel('Temperature (deg C)')
Explanation: We can, for example, plot the current temperature versus latitude:
End of explanation
#  The dictionary of sub-processes:
mymodel.subprocess
Explanation: It is based on a very general concept of a model as a collection of individual, interacting processes. 
Every model in climlab is a collection of individual, interacting physical processes. 
In this case, the EBM_annual object has a list of 4 sub-processes. Each one is basically one of the terms in our energy budget equation above:
diffusion
Longwave radiation ($A+BT$)
Insolation (annual mean)
Albedo (in this case, a fixed spatial variation)
To see this explicitly:
End of explanation
#  Make a copy of the original temperature array
initial = mymodel.Ts.copy()
#  Take a single timestep forward!
mymodel.step_forward()
#  Check out the difference
print mymodel.Ts - initial
Explanation: So what does it do?
Well for one thing, we can step the model forward in time. This is basically just like we did manually with our zero-dimensional model
End of explanation
#  How long is a single timestep?
mymodel.timestep
Explanation: Looks like the temperature got a bit colder near the equator and warmer near the poles
End of explanation
mymodel.param
Explanation: This value is in seconds. It is actually 1/90th of a year (so, to step forward one year, we need 90 individual steps). This is a default value -- we could change it if we wanted to.
Actually we can see a bunch of important parameters for this model all together in a dictionary called param:
End of explanation
#  Each climlab model has a dictionary called diagnostics.
#  Let's look at the names of the fields in this dicionary:
mymodel.diagnostics.keys()
#  We can access individual fields either through standard dictionary notation:
mymodel.diagnostics['ASR']
#  Or using the more interactive 'dot' notation:
mymodel.ASR
#  Let's use the diagnostics to make a plot of the current state of the radiation
plt.plot(mymodel.lat, mymodel.ASR, label='ASR')
plt.plot(mymodel.lat, mymodel.OLR, label='OLR')
plt.xlabel('Latitude')
plt.ylabel('W/m2')
plt.legend()
plt.grid()
Explanation: Accessing the model diagnostics
Each model can have lots of different diagnostic quantities. This means anything that can be calculated from the current model state -- in our case, from the temperature.
For the EBM, this includes (among other things) the OLR and the ASR.
End of explanation
climlab.global_mean(mymodel.net_radiation)
climlab.global_mean(mymodel.Ts)
Explanation: This plot shows that $ASR > OLR$ (system is gaining extra energy) across the tropics, and $ASR < OLR$ (system is losing energy) near the poles.
It's interesting to ask how close this model is to global energy balance.
For that, we have to take an area-weighted global average.
End of explanation
#  Loop 90 times for 1 year of simulation
for n in range(90):
    mymodel.step_forward()
print 'Global mean temperature is %0.1f degrees C.' %climlab.global_mean(mymodel.Ts)
print 'Global mean energy imbalance is %0.1f W/m2.' %climlab.global_mean(mymodel.net_radiation)
Explanation: Running the model out to equilibrium
We write a loop to integrate the model through many timesteps, and then check again for energy balance:
End of explanation
mymodel.integrate_years(1.)
print 'Global mean temperature is %0.1f degrees C.' %climlab.global_mean(mymodel.Ts)
print 'Global mean energy imbalance is %0.1f W/m2.' %climlab.global_mean(mymodel.net_radiation)
Explanation: Since there is still a significant imbalance, we are not yet at equilibrium. We should step forward again.
This time, let's use a convenient built-in function instead of writing our own loop:
End of explanation
plt.plot(mymodel.lat, mymodel.Ts)
plt.xlabel('Latitude')
plt.ylabel('Temperature (deg C)')
plt.grid()
plt.plot(mymodel.lat_bounds, mymodel.heat_transport())
plt.xlabel('Latitude')
plt.ylabel('PW')
plt.grid()
Explanation: We are now quite close to equilibrium. Let's make some plots
End of explanation
def ebm_plot(e, return_fig=False):    
    templimits = -20,32
    radlimits = -340, 340
    htlimits = -6,6
    latlimits = -90,90
    lat_ticks = np.arange(-90,90,30)
    
    fig = plt.figure(figsize=(8,12))
    ax1 = fig.add_subplot(3,1,1)
    ax1.plot(e.lat, e.Ts)
    ax1.set_ylim(templimits)
    ax1.set_ylabel('Temperature (deg C)')
    
    ax2 = fig.add_subplot(3,1,2)
    ax2.plot(e.lat, e.ASR, 'k--', label='SW' )
    ax2.plot(e.lat, -e.OLR, 'r--', label='LW' )
    ax2.plot(e.lat, e.net_radiation, 'c-', label='net rad' )
    ax2.plot(e.lat, e.heat_transport_convergence(), 'g--', label='dyn' )
    ax2.plot(e.lat, e.net_radiation.squeeze() + e.heat_transport_convergence(), 'b-', label='total' )
    ax2.set_ylim(radlimits)
    ax2.set_ylabel('Energy budget (W m$^{-2}$)')
    ax2.legend()
    
    ax3 = fig.add_subplot(3,1,3)
    ax3.plot(e.lat_bounds, e.heat_transport() )
    ax3.set_ylim(htlimits)
    ax3.set_ylabel('Heat transport (PW)')
    
    for ax in [ax1, ax2, ax3]:
        ax.set_xlabel('Latitude')
        ax.set_xlim(latlimits)
        ax.set_xticks(lat_ticks)
        ax.grid()
    
    if return_fig:
        return fig
ebm_plot(mymodel)
Explanation: Let's make some nice plots of all the terms in the energy budget.
We'll define a reusable function to make the plots.
End of explanation
model2 = climlab.EBM_annual(num_lat=30)
print model2
#  The default initial temperature
print model2.Ts
#  Now let's change it to be 15 degrees everywhere
model2.Ts[:] = 15.
print model2.Ts
model2.compute_diagnostics()
ebm_plot(model2)
Explanation: What if we start from a very different initial temperature?
End of explanation
#  Code to generate individual frames for the animation
#  You need to specify a valid path for your computer, or else this won't work
#  Uncomment the code below to re-create the frames from the animation
#nsteps = 90
#for n in range(nsteps):
#    fig = ebm_plot(model2, return_fig=True)
#    
#    filepath = '/Users/br546577/Desktop/temp/ebm_animation_frames/'
#    filename = filepath + 'ebm_animation' + str(n).zfill(4) + '.png'
#    fig.savefig(filename)
#    plt.close()
#    
#    model2.step_forward()
Explanation: Why is the heat transport zero everywhere?
Now let's look at an animation showing how all the terms in the energy budget evolve from this initial condition toward equilibrium.
The code below generates a series of snapshots, that we can stitch together into an animation with a graphics program.
End of explanation
<END_TASK>