In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
Open the output from our control simulation with the slab ocean version of the CESM:
In [2]:
## To read data over the internet
control_filename = 'som_1850_f19.cam.h0.clim.nc'
datapath = 'http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/Brian+Rose/CESM+runs/'
endstr = '/entry.das'
control = nc.Dataset( datapath + 'som_1850_f19/' + control_filename + endstr )
## To read from a local copy of the file
## (just a small subset of the total list of variables, to save disk space)
#ontrol_filename = 'som_1850_f19.cam.h0.clim_subset.nc'
#control = nc.Dataset( control_filename )
The full file from the online server contains many many variables, describing all aspects of the model climatology.
Whether we see a long list or a short list in the following code block depends on whether we are reading the full output file or the much smaller subset:
In [3]:
for v in control.variables: print v
Today we need just a few of these variables:
TS: the surface temperatureFLNT: the longwave radiation at the top of the atmosphere (i.e. what we call the OLR)FSNT: the net shortwave radiation at the top of the atmosphere (i.e. what we call the ASR)FLNTC: the clear-sky OLRFSNTC: the clear-sky ASRTake a look at some of the meta-data for these fields:
In [4]:
for field in ['TS', 'FLNT', 'FSNT', 'FLNTC', 'FSNTC']:
print control.variables[field]
Each one of these variables has dimensions (12, 96, 144), which corresponds to time (12 months), latitude and longitude.
Take a look at one of the coordinate variables:
In [5]:
print control.variables['lat']
Now let's load in the coordinate data, to use later for plotting:
In [6]:
lat = control.variables['lat'][:]
lon = control.variables['lon'][:]
print lat
In [7]:
# A re-usable function to make a map of a 2d field on a latitude / longitude grid
def make_map(field_2d):
# Make a filled contour plot
fig = plt.figure(figsize=(10,5))
cax = plt.contourf(lon, lat, field_2d)
# draw a single contour to outline the continents
plt.contour( lon, lat, control.variables['LANDFRAC'][0,:,:], [0.5], colors='k')
plt.xlabel('Longitude (degrees east)')
plt.ylabel('Latitude (degrees north)')
plt.colorbar(cax)
In [8]:
# Here is a convenient function that takes the name of a variable in our CESM output
# and make a map of its annual average
def map_this(fieldname, dataset=control):
field = dataset.variables[fieldname][:]
field_annual = np.mean(field, axis=0)
make_map(field_annual)
In [9]:
# Use this function to make a quick map of the annual average surface temperature:
map_this('TS')
In [10]:
# The lat/lon dimensions after taking the time average:
TS_annual = np.mean(control.variables['TS'][:], axis=0)
TS_annual.shape
Out[10]:
Define a little re-usable function to take the global average of any of these fields:
In [11]:
def global_mean(field_2d):
'''This function takes a 2D array on a regular latitude-longitude grid
and returns the global area-weighted average'''
zonal_mean = np.mean(field_2d, axis=1)
return np.average(zonal_mean, weights=np.cos(np.deg2rad(lat)))
In [12]:
# Again, a convenience function that takes just the name of the model output field
# and returns its time and global average
def global_mean_this(fieldname, dataset=control):
field = dataset.variables[fieldname][:]
field_annual = np.mean(field, axis=0)
return global_mean(field_annual)
Now compute the global average surface temperature in the simulation:
In [13]:
global_mean_this('TS')
Out[13]:
The model simulates cloud amount in every grid box. The cloud field is thus 4-dimensional:
In [14]:
# This field is not included in the small subset file
# so this will only work if you are reading the full file from the online server
control.variables['CLOUD']
Out[14]:
To simplify things we can just look at the total cloud cover, integrated from the surface to the top of the atmosphere:
In [15]:
control.variables['CLDTOT']
Out[15]:
In [16]:
map_this('CLDTOT')
Which parts of Earth are cloudy and which are not? (at least in this simulation)
In [17]:
# To get you started, here is the ASR
map_this('FSNT')
In [18]:
map_this('FLNT')
In [19]:
net_radiation = np.mean(control.variables['FSNT'][:] - control.variables['FLNT'][:], axis=0)
make_map(net_radiation)
In [20]:
global_mean(net_radiation)
Out[20]:
In [ ]:
In [ ]:
In [ ]:
How much CO2 was in the atmosphere for the control simulation?
This information is available in the full output file (this won't work with the local subset file):
In [ ]:
# The meta-data:
control.variables['co2vmr']
In [ ]:
# The data themselves, expressed in ppm:
control.variables['co2vmr'][:] * 1E6
Answer: the CO2 concentration is 284.7 ppm in the control simulation.
Now we want to see how the climate changes in the CESM when we double CO2 and run it out to equilibrium.
I have done this. Because we are using a slab ocean model, it reaches equilibrium after just a few decades.
Let's now open up the output file from the 2xCO2 scenario:
In [ ]:
## To read data over the internet
# doubleCO2_filename = 'som_1850_2xCO2.cam.h0.clim.nc'
# doubleCO2 = nc.Dataset( datapath + 'som_1850_f19/' + doubleCO2_filename + endstr )
## To read from a local copy of the file
## (just a small subset of the total list of variables, to save disk space)
doubleCO2_filename = 'som_1850_2xCO2.cam.h0.clim_subset.nc'
doubleCO2 = nc.Dataset( doubleCO2_filename )
This file has all the same fields as control, but they reflect the new equilibrium climate after doubling CO2.
Let's verify the CO2 amount:
In [ ]:
doubleCO2.variables['co2vmr'][:] * 1E6
So the CO2 concentration is now 569.4 ppm.
In [ ]:
In [ ]:
Are the clouds warming or cooling the climate?
How has CRE changed compared to the control climate? (i.e. is the net effect larger or smaller than it was before)
Can you infer whether the cloud feedback is positive or negative in the CESM?
In [ ]:
In [ ]:
In [ ]: