Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

... work in progress.

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs

The CFSR climatology datasets

We pre-computed 30-year seasonal and monthly climatologies from the 6-hourly CFSR dataset, see Computing seasonal and monthly means from CFSR data.

Here we’ll use those datasets to make some nice reference figures.

local_path = '/nfs/roselab_rit/data/cfsr_climatology/'

path = local_path

Load the seasonal data:

cfsr_seas = xr.open_mfdataset(path + '*' + '.seas_clim.0p5.nc')
cfsr_seas
Loading...

and the monthly data:

cfsr_mon = xr.open_mfdataset(path + '*' + '.mon_clim.0p5.nc')
cfsr_mon
Loading...

Sea level pressure

For many fields, we will make three individual plots:

  • annual mean

  • DJF

  • JJA

Here’s an example for the sea level pressure. I’m choosing a particular map projection here. There are lots of other choices, see the Cartopy documentation.

lon = cfsr_seas.lon
lat = cfsr_seas.lat

levels = np.arange(982., 1040., 4.)

fig, axes = plt.subplots(3, figsize=(10,15), subplot_kw={'projection': ccrs.Robinson()})

for ax in axes:
    ax.coastlines()
    
ax = axes[0]
CS = ax.contour(lon, lat, cfsr_seas.pmsl.mean(dim='season') / 100., levels=levels, linewidths=0.8, 
               transform=ccrs.PlateCarree())
ax.clabel(CS, CS.levels, inline=True, fontsize=8)
ax.set_title('Annual mean')

ax = axes[1]
CS = ax.contour(lon, lat, cfsr_seas.pmsl.sel(season='DJF') / 100., levels=levels, linewidths=0.8, 
               transform=ccrs.PlateCarree())
ax.clabel(CS, CS.levels, inline=True, fontsize=8)
ax.set_title('DJF')

ax = axes[2]
CS = ax.contour(lon,lat, cfsr_seas.pmsl.sel(season='JJA') / 100., levels=levels, linewidths=0.8, 
               transform=ccrs.PlateCarree())
ax.clabel(CS, CS.levels, inline=True, fontsize=8)
ax.set_title('JJA')
<Figure size 1000x1500 with 3 Axes>

List of figures

The following figures appear in the slide deck we saw in class, which you can find at this link (courtesy of Professor Brian Tang, University at Albany).

Our objective here is to build recipes to make all (or at least most of) these figures from the CFSR climatology.

  1. SLP this one is already done

    • annual mean

    • DJF

    • JJA

  2. 1000 hPa streamlines

    • annual mean

    • DJF

    • JJA

  3. 1000 hPa zonal wind

    • annual mean

    • DJF

    • JJA

  4. 1000 hPa meridional wind

    • annual mean

    • DJF

    • JJA

  5. 1000 hPa temperature

    • annual mean

  6. 1000 hPa zonal temperature anomaly

    • DJF

    • JJA

  7. 200 hPa streamlines

    • annual mean

    • DJF

    • JJA

  8. 200 hPa zonal wind

    • annual mean

    • DJF

    • JJA

  9. 200 hPa meridional wind

    • annual mean

    • DJF

    • JJA

  10. 500 hPa vertical velocity (omega)

    • annual mean

    • DJF

    • JJA

  11. 500 hPa relative humidity

    • annual mean

    • DJF

    • JJA

  12. 850 hPa moist static energy

    • annual mean

  13. 850 hPa moist static energy zonal anomaly

    • DJF

    • JJA

  14. Surface enthalpy flux

    • annual mean

    • DJF

    • JJA

  15. Schematic: Lagrangian vs Eulerian overturning

Mass overturning circulation

Basic ideas and theory

An element of northward mass flux (in units of kg s−1{^-1}) across latitude line ϕ\phi is

dψ=acos⁡ϕdλ v dpgd\psi = a \cos\phi d\lambda ~ v ~ \frac{dp}{g}
(1)

So the total northward mass flux across ϕ\phi above pressure level pp is given by

ψ(p)=∫02π∫0pacos⁡ϕdλ v dpg=2πacos⁡ϕg∫0p[v]dp\begin{align*} \psi(p) &= \int_0^{2\pi} \int_0^p a \cos\phi d\lambda ~ v ~ \frac{dp}{g} \\ &= 2\pi a \frac{\cos\phi}{g} \int_0^p [v] dp \end{align*}
(2)

The zonally- and time-averaged mass conservation in pressure coordinates is

1acos⁡ϕ∂∂ϕ([v‾]cos⁡ϕ)+∂[ω‾]∂p=0\frac{1}{a \cos\phi} \frac{\partial}{\partial \phi} \left( [\overline{v}] \cos\phi \right) + \frac{\partial [\overline{\omega}]}{\partial p} = 0
(3)

In other words, the vector with components [v‾],[ω‾][\overline{v}], [\overline{\omega}] is non-divergent. So we can represent it with a scalar streamfunction ψ\psi, the mass transport streamfunction:

2πacos⁡ϕ[v‾]=∂ψ∂p2πa2cos⁡ϕ([ω‾])=−∂ψ∂ϕ\begin{align*} 2\pi a \cos\phi [\overline{v}] &= \frac{\partial \psi}{\partial p} \\ 2\pi a^2 \cos\phi \left( [\overline{\omega}]\right) = -\frac{\partial \psi}{\partial \phi} \end{align*}
(4)

consistent with the definition of ψ\psi above!

It it thus possible to compute and plot the mass streamfunction from just the zonal wind data, using

ψ=2πacos⁡ϕg∫0p[v‾]dp\psi = \frac{2\pi a \cos\phi}{g} \int_0^p [\overline{v}] dp
(5)

Contours of ψ\psi can then be interpreted as streamlines for the mean meridional circulation, that is, the zonal- and time-averaged overturning of mass in the latitude-pressure plane.

Computing the overturning streamfunction from meridonal velocities

We’ll start by taking averages of the vv data:

vbracket = cfsr_seas.v.mean(dim='lon')
vbracket
Loading...

It turns out that the top-most pressure level is missing in the data, e.g.:

vbracket.isel(lev=0).sel(season='DJF').values
array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan], dtype=float32)

We can’t integrate over nans so we’ll replace those with zeros instead -- which will make no contribution to the integral.

vbracket_zeros = vbracket.fillna(0.)

Now let’s compute the integral!

(I’m calling out to an external package to get values of the constants aa (Earth radius) and gg (acceleration due to gravity))

coslat = np.cos(np.deg2rad(cfsr_seas.lat))
from climlab.utils.constants import a, g

# Use units of 10^9 kg/s
psi = 2*np.pi*a*coslat/g * vbracket_zeros.cumulative_integrate(coord='lev') / 1E9 * 100
# Need to multiply by 100 because the 'lev' coordinate is in units of hPa, not Pa
psi
Loading...

Now we’ve got the streamfunction defined for four seasons. Time to make plots!

Plotting the mass overturning

psi_levels = np.arange(-230, 240, 10.)

fig, axes = plt.subplots(3, figsize=(8,10))

CS = psi.mean(dim='season').plot.contour(ax=axes[0],
                                         x='lat', 
                                         yincrease=False,
                                         levels=psi_levels,
                                        )
axes[0].clabel(CS, CS.levels, inline=True, fontsize=8)
axes[0].set_title('Annual mean mass streamfunction (10$^9$ kg s$^{-1}$)')

CS = psi.sel(season='DJF').plot.contour(ax=axes[1],
                                         x='lat', 
                                         yincrease=False,
                                         levels=psi_levels,
                                        )
axes[1].clabel(CS, CS.levels, inline=True, fontsize=8)
axes[1].set_title('DJF mass streamfunction (10$^9$ kg s$^{-1}$)')

CS = psi.sel(season='JJA').plot.contour(ax=axes[2],
                                         x='lat', 
                                         yincrease=False,
                                         levels=psi_levels,
                                        )
axes[2].clabel(CS, CS.levels, inline=True, fontsize=8)
axes[2].set_title('JJA mass streamfunction (10$^9$ kg s$^{-1}$)')
<Figure size 800x1000 with 3 Axes>
Lectures
Computing seasonal and monthly means from CFSR data
Lectures
Decompositions of the circulation