3. Analyze CESM Output

paleoCAMP logo

3. Analyze CESM Output#

Tutorials at the 2026 paleoCAMP | June 15–June 29, 2026

Jiang Zhu
jiangzhu@ucar.edu
Climate & Global Dynamics Laboratory
NSF National Center for Atmospheric Research

Learning Objectives:#

  • Learn to use the NCAR JupyterHub for data access and analysis

  • Learn to read and examine NetCDF files using Xarray

  • Learn techniques to make basic visualization of CESM output of temperature, precipitation, and sea-surface temperature

Time to learn: 60 minutes

How to get started?

  1. Launch a JupyterHub server: https://jupyterhub.hpc.ucar.edu/

  • Sign in with your username and password,passcode

    • password is the CIT password you set up with CISL

    • passcode is the six digits on your DUO app

    • NOTE the comma , between password and passcode

DUO Passcode

Figure: Screenshot of DUO app

  1. Start the JupyterHub. From the dropdown menu, select

  • Casper PBS Batch for Resource Selection

  • casper for Queue or Reservation

  • UAZN0052 for Project Account

  • 06:00:00 for Wall Time

Jupyter Login

Figure: Jupyterhub Casper

  1. Launch a terminal: click File, New, and Terminal

Launch Terminal

Figure: Launch a Terminal

  1. Type git clone https://github.com/jiang-zhu/paleocamp2026.git

  2. Find and open 3_analyze_CESM_output.ipynb from the left sidebar

  3. Select a kernel from top right corner NPL 2025a or NPL 2026a Select kernel

Figure: Select an NPL kernel

Load Python packages

import os
import glob
from datetime import timedelta

import xarray as xr
import numpy as np
import pandas as pd

import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point

# geocat is used for interpolation of the atmosphere output
from geocat.comp import interpolation

# xesmf is used for regridding the ocean output
import xesmf

# hvplot provides interactive plots
import hvplot.xarray

import warnings
warnings.simplefilter("ignore", FutureWarning)

Analysis 1: plot solar insolation in the MH to further validate the simulation#

Load data#

  • I have both the piControl and midHolocene simulations finished.

  • I see 12-month history files for each simulation in my Archive Directory.

!ls /glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/*cam.h0.0001*
!ls /glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/*cam.h0.0001*

/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-01.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-02.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-03.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-04.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-05.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-06.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-07.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-08.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-09.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-10.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-11.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-12.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-01.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-02.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-03.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-04.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-05.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-06.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-07.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-08.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-09.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-10.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-11.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-12.nc

  • Use glob and os.path.join to obtain the file names including path.

  • Use wildcard * to catch all files in the atmosphere history directory (see Wildcard in 0.1_demo_unix.ipynb).

storage_dir = '/glade/derecho/scratch/jiangzhu/archive/'
hist_dir = 'atm/hist'

case_PI = 'b.e21.B1850.f19_g17.piControl.001'
case_MH = 'b.e21.B1850.f19_g17.midHolocene.001'

files_PI = glob.glob(os.path.join(storage_dir, case_PI, hist_dir, '*.cam.h0.0001*'))
files_MH = glob.glob(os.path.join(storage_dir, case_MH, hist_dir, '*.cam.h0.0001*'))

print('List of files for PI')
print(*files_PI, sep='\n')

print('List of files for MH')
print(*files_MH, sep='\n')

List of files for PI
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-10.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-02.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-01.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-11.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-06.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-03.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-05.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-12.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-04.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-09.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-07.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.piControl.001/atm/hist/b.e21.B1850.f19_g17.piControl.001.cam.h0.0001-08.nc
List of files for MH
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-04.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-06.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-07.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-03.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-11.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-08.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-01.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-05.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-02.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-09.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-12.nc
/glade/derecho/scratch/jiangzhu/archive/b.e21.B1850.f19_g17.midHolocene.001/atm/hist/b.e21.B1850.f19_g17.midHolocene.001.cam.h0.0001-10.nc

  • Use xr.open_mfdataset to open multiple files at once in parallel (12 months in this case)

  • We will keep using the Xarray Datasets ds_MH and ds_PI throughout the Notebook

ds_PI = xr.open_mfdataset(
    files_PI,
    combine="by_coords",
    data_vars="minimal",
    coords="minimal",
    compat="override",
)

ds_MH = xr.open_mfdataset(
    files_MH,
    combine="by_coords",
    data_vars="minimal",
    coords="minimal",
    compat="override",
)

ds_PI

<xarray.Dataset> Size: 541MB
Dimensions:       (zlon: 1, nbnd: 2, lat: 96, lev: 26, ilev: 27, time: 12,
                   lon: 144)
Coordinates:
  * zlon          (zlon) float64 8B 0.0
  * lat           (lat) float64 768B -90.0 -88.11 -86.21 ... 86.21 88.11 90.0
  * lev           (lev) float64 208B 3.545 7.389 13.97 ... 929.6 970.6 992.6
  * ilev          (ilev) float64 216B 2.194 4.895 9.882 ... 956.0 985.1 1e+03
  * time          (time) object 96B 0001-02-01 00:00:00 ... 0002-01-01 00:00:00
  * lon           (lon) float64 1kB 0.0 2.5 5.0 7.5 ... 350.0 352.5 355.0 357.5
Dimensions without coordinates: nbnd
Data variables: (12/122)
    zlon_bnds     (zlon, nbnd) float64 16B dask.array<chunksize=(1, 2), meta=np.ndarray>
    gw            (lat) float64 768B dask.array<chunksize=(96,), meta=np.ndarray>
    hyam          (lev) float64 208B dask.array<chunksize=(26,), meta=np.ndarray>
    hybm          (lev) float64 208B dask.array<chunksize=(26,), meta=np.ndarray>
    P0            float64 8B ...
    hyai          (ilev) float64 216B dask.array<chunksize=(27,), meta=np.ndarray>
    ...            ...
    VU            (time, lev, lat, lon) float32 17MB dask.array<chunksize=(1, 26, 96, 144), meta=np.ndarray>
    VV            (time, lev, lat, lon) float32 17MB dask.array<chunksize=(1, 26, 96, 144), meta=np.ndarray>
    Vzm           (time, ilev, lat, zlon) float32 124kB dask.array<chunksize=(1, 27, 96, 1), meta=np.ndarray>
    WTHzm         (time, ilev, lat, zlon) float32 124kB dask.array<chunksize=(1, 27, 96, 1), meta=np.ndarray>
    Wzm           (time, ilev, lat, zlon) float32 124kB dask.array<chunksize=(1, 27, 96, 1), meta=np.ndarray>
    Z3            (time, lev, lat, lon) float32 17MB dask.array<chunksize=(1, 26, 96, 144), meta=np.ndarray>
Attributes:
    Conventions:       CF-1.0
    source:            CAM
    case:              b.e21.B1850.f19_g17.piControl.001
    logname:           jiangzhu
    host:              derecho3
    initial_file:      /glade/campaign/cesm/cesmdata/inputdata/atm/cam/inic/f...
    topography_file:   /glade/campaign/cesm/cesmdata/inputdata/atm/cam/topo/f...
    model_doi_url:     https://doi.org/10.5065/D67H1H0V
    time_period_freq:  month_1

  • Click time in Coordinates of the above data array. It starts with 0001-02-01, which is wrong

  • We need to fix this CESM issue to get the correct time, i.e., January to December of year 1

ds_PI['time'] = ds_PI.time.get_index('time') - timedelta(days=15)
ds_MH['time'] = ds_MH.time.get_index('time') - timedelta(days=15)

ds_PI

# Explore the dataset and find the Solar insolation, SOLIN

<xarray.Dataset> Size: 541MB
Dimensions:       (zlon: 1, nbnd: 2, lat: 96, lev: 26, ilev: 27, time: 12,
                   lon: 144)
Coordinates:
  * zlon          (zlon) float64 8B 0.0
  * lat           (lat) float64 768B -90.0 -88.11 -86.21 ... 86.21 88.11 90.0
  * lev           (lev) float64 208B 3.545 7.389 13.97 ... 929.6 970.6 992.6
  * ilev          (ilev) float64 216B 2.194 4.895 9.882 ... 956.0 985.1 1e+03
  * time          (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
  * lon           (lon) float64 1kB 0.0 2.5 5.0 7.5 ... 350.0 352.5 355.0 357.5
Dimensions without coordinates: nbnd
Data variables: (12/122)
    zlon_bnds     (zlon, nbnd) float64 16B dask.array<chunksize=(1, 2), meta=np.ndarray>
    gw            (lat) float64 768B dask.array<chunksize=(96,), meta=np.ndarray>
    hyam          (lev) float64 208B dask.array<chunksize=(26,), meta=np.ndarray>
    hybm          (lev) float64 208B dask.array<chunksize=(26,), meta=np.ndarray>
    P0            float64 8B ...
    hyai          (ilev) float64 216B dask.array<chunksize=(27,), meta=np.ndarray>
    ...            ...
    VU            (time, lev, lat, lon) float32 17MB dask.array<chunksize=(1, 26, 96, 144), meta=np.ndarray>
    VV            (time, lev, lat, lon) float32 17MB dask.array<chunksize=(1, 26, 96, 144), meta=np.ndarray>
    Vzm           (time, ilev, lat, zlon) float32 124kB dask.array<chunksize=(1, 27, 96, 1), meta=np.ndarray>
    WTHzm         (time, ilev, lat, zlon) float32 124kB dask.array<chunksize=(1, 27, 96, 1), meta=np.ndarray>
    Wzm           (time, ilev, lat, zlon) float32 124kB dask.array<chunksize=(1, 27, 96, 1), meta=np.ndarray>
    Z3            (time, lev, lat, lon) float32 17MB dask.array<chunksize=(1, 26, 96, 144), meta=np.ndarray>
Attributes:
    Conventions:       CF-1.0
    source:            CAM
    case:              b.e21.B1850.f19_g17.piControl.001
    logname:           jiangzhu
    host:              derecho3
    initial_file:      /glade/campaign/cesm/cesmdata/inputdata/atm/cam/inic/f...
    topography_file:   /glade/campaign/cesm/cesmdata/inputdata/atm/cam/topo/f...
    model_doi_url:     https://doi.org/10.5065/D67H1H0V
    time_period_freq:  month_1

Compute the zonal mean of solar insolation and make plots#

  • Examine a variable including its dimension, long name, and units

  • Use the .isel method to select a single month or multiple months

  • Use the Xarray plotting functionality to make a simple plot

  • Use the hvplot to make an interactive plot

  • Use .mean('lon') to get the zonal mean

ds_PI.SOLIN

<xarray.DataArray 'SOLIN' (time: 12, lat: 96, lon: 144)> Size: 664kB
dask.array<concatenate, shape=(12, 96, 144), dtype=float32, chunksize=(1, 96, 144), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    Sampling_Sequence:  rad_lwsw
    units:              W/m2
    long_name:          Solar insolation
    cell_methods:       time: mean
ds_PI.SOLIN.isel(time=5).plot(size=2)

<matplotlib.collections.QuadMesh at 0x1531d9cfc830>
_images/f121413128e41f4043b165f5608fbf04255489d316f3d61ba958f5c535844b6f.png 
# Easier to find value for SNARL (37.6N, -118.8E)

ds_PI.SOLIN.isel(time=5).hvplot(coastline=True, cmap='viridis')

ds_PI.SOLIN

<xarray.DataArray 'SOLIN' (time: 12, lat: 96, lon: 144)> Size: 664kB
dask.array<concatenate, shape=(12, 96, 144), dtype=float32, chunksize=(1, 96, 144), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    Sampling_Sequence:  rad_lwsw
    units:              W/m2
    long_name:          Solar insolation
    cell_methods:       time: mean
# Compute zonal mean and make plot of piControl
ds_PI.SOLIN.mean('lon').plot.contourf(x='time', y='lat', figsize=(4, 2))
plt.xlabel('time (year-month)')

Text(0.5, 0, 'time (year-month)')
_images/247419dbc81bd2b2b3ca6a94e9a063e7367b88a3b1a05d06b4cfe5d0508c4a81.png 
# Add your code to compute zonal mean and make plot of midHolocene

# Add your code to compute and plot the difference: midHolocen - piControl zonal mean
Click here for the solution
Copy and paste the code into the above cell 
ds_MH.SOLIN.mean('lon').plot.contourf(x='time', y='lat', figsize=(4, 2))

(ds_MH.SOLIN - ds_PI.SOLIN).mean('lon').plot.contourf(
    x='time', y='lat', figsize=(4, 2), levels=np.linspace(-30, 30, 21))

Small group discussion#

  • Which orbital parameters are different at the middle Holocene (6ka BP)?

  • How does the orbital parameter impact the top-of-atmosphere shortwave radiation (solar insolation)?

  • Do the results look correct? You can compare your results with Figure 3b of Otto-Bliesner et al., (2017)

Analysis 2: Does the Earth receive more radiation during the mid-Holocene?#

  • We compute the global annual mean solar insolation to answer this question.

  • Importantly, we need to weight grid cells by their area (equivalent to cosine of latitude), using the .weighted method

Note that lat is in degrees

ds_PI.lat

<xarray.DataArray 'lat' (lat: 96)> Size: 768B
array([-90.      , -88.105263, -86.210526, -84.315789, -82.421053, -80.526316,
       -78.631579, -76.736842, -74.842105, -72.947368, -71.052632, -69.157895,
       -67.263158, -65.368421, -63.473684, -61.578947, -59.684211, -57.789474,
       -55.894737, -54.      , -52.105263, -50.210526, -48.315789, -46.421053,
       -44.526316, -42.631579, -40.736842, -38.842105, -36.947368, -35.052632,
       -33.157895, -31.263158, -29.368421, -27.473684, -25.578947, -23.684211,
       -21.789474, -19.894737, -18.      , -16.105263, -14.210526, -12.315789,
       -10.421053,  -8.526316,  -6.631579,  -4.736842,  -2.842105,  -0.947368,
         0.947368,   2.842105,   4.736842,   6.631579,   8.526316,  10.421053,
        12.315789,  14.210526,  16.105263,  18.      ,  19.894737,  21.789474,
        23.684211,  25.578947,  27.473684,  29.368421,  31.263158,  33.157895,
        35.052632,  36.947368,  38.842105,  40.736842,  42.631579,  44.526316,
        46.421053,  48.315789,  50.210526,  52.105263,  54.      ,  55.894737,
        57.789474,  59.684211,  61.578947,  63.473684,  65.368421,  67.263158,
        69.157895,  71.052632,  72.947368,  74.842105,  76.736842,  78.631579,
        80.526316,  82.421053,  84.315789,  86.210526,  88.105263,  90.      ])
Coordinates:
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
Attributes:
    long_name:  latitude
    units:      degrees_north

Example calculation for the piControl#

coslat = np.cos(np.deg2rad(ds_PI.lat))

SOLIN_PI = ds_PI.SOLIN.weighted(coslat).mean(('lat', 'lon', 'time'))

print("The global annual mean insolation of PI: ", SOLIN_PI.values, "W/m2")

The global annual mean insolation of PI:  340.32719678192876 W/m2

If we don’t apply the area weights, we will get a wrong answer#

SOLIN_PI_unweighted = ds_PI.SOLIN.mean(('lat', 'lon', 'time'))

print("Only if we forgot to properly weight the values by area (cosine(latitude)):",
      SOLIN_PI_unweighted.values, "W/m2")

Only if we forgot to properly weight the values by area (cosine(latitude)): 297.03168 W/m2

Add your own calculation for the midHolocene and compute the difference between midHolocene and piControl#

  • SOLIN_MH = ?

  • SOLIN_MH - SOLIN_PI = ?

Click here for the solution Copy and paste the code into the above cell 
coslat = np.cos(np.deg2rad(ds_MH.lat))

SOLIN_MH = ds_MH.SOLIN.weighted(coslat).mean(('lat', 'lon', 'time'))
print("The global annual mean insolation of MH: ", SOLIN_MH.values, "W/m2")
print("Difference, MH - PI: ", SOLIN_MH.values-SOLIN_PI.values, "W/m2")

Small group discussion#

  • Jess taught us in her lecture that obliquity and precession do not impact the global annual mean solar insolation, but why do we see different value in our midHolonece?

Click here for the solution

  • The code below should produce a larger value for midHolocene by ~0.012 W m-2 (0.004%).

  • The larger solar insolation is very likely due to the higher eccentricity in midHolocene (0.018682 vs 0.01676429)!

coslat = np.cos(np.deg2rad(ds_MH.lat))
month_length = ds_MH.time.dt.days_in_month

SOLIN_PI = ds_PI.SOLIN.weighted(month_length * coslat).mean(('time', 'lat', 'lon'))
SOLIN_MH = ds_MH.SOLIN.weighted(month_length * coslat).mean(('time', 'lat', 'lon'))

print("The global annual mean insolation of PI: ", SOLIN_PI.values, "W/m2")
print("The global annual mean insolation of MH: ", SOLIN_MH.values, "W/m2")
print("Difference, MH - PI: ", SOLIN_MH.values - SOLIN_PI.values, "W/m2")

Small group discussion#

  • Should the mid-Holocene be warmer or colder than the preindustrial (given that difference in solar insolation is only 0.012 W m-2, 0.004%)?

  • Holocene Temperature Conundrum: geological records suggest that mid-Holocene was warmer by ~0.5℃ than the preindustrial, while climate models suggest a colder mid-Holocene by ~0.5℃.

Analysis 3: how does orbital forcing impact the ITCZ and monsoon precipitation during the midHolocene?#

  • In CESM, total precipitation is computed as prec = PRECC + PRECL (sum of convective and large-scale precipitation; recall that climate models have to parameterize convection!)

  • We convert the units from m/s into mm/day

  • We use .isel(time=slice(5, 8)) to select the June, July, and August values (recall that Python uses 0-based ordering)

  • Let’s use Cartopy and Matplotlib with a Robinson projection (instead of using the simple Xarray.plot())

  • We use add_cyclic_point to get rid of the “white strip” in the plot

ds_PI.PRECC

<xarray.DataArray 'PRECC' (time: 12, lat: 96, lon: 144)> Size: 664kB
dask.array<concatenate, shape=(12, 96, 144), dtype=float32, chunksize=(1, 96, 144), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    units:         m/s
    long_name:     Convective precipitation rate (liq + ice)
    cell_methods:  time: mean
ds_PI.PRECL

<xarray.DataArray 'PRECL' (time: 12, lat: 96, lon: 144)> Size: 664kB
dask.array<concatenate, shape=(12, 96, 144), dtype=float32, chunksize=(1, 96, 144), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    units:         m/s
    long_name:     Large-scale (stable) precipitation rate (liq + ice)
    cell_methods:  time: mean
m_p_s_to_mm_p_day = 86400 * 1000

prec_PI = (ds_PI.PRECC + ds_PI.PRECL).isel(time=slice(5, 8)).mean('time') * m_p_s_to_mm_p_day

prec_MH = (ds_MH.PRECC + ds_MH.PRECL).isel(time=slice(5, 8)).mean('time') * m_p_s_to_mm_p_day

prec_PI

<xarray.DataArray (lat: 96, lon: 144)> Size: 55kB
dask.array<mul, shape=(96, 144), dtype=float32, chunksize=(96, 144), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    units:         m/s
    cell_methods:  time: mean

Plot total precipitation of piControl#

fig, ax = plt.subplots(figsize=(3, 1.5), subplot_kw={
    'projection': ccrs.Robinson(central_longitude=210)})

p1 = ax.contourf(prec_PI.lon, prec_PI.lat, prec_PI,
                 cmap='YlGnBu',
                 levels=np.linspace(0, 17, 18),
                 extend='both',
                 transform=ccrs.PlateCarree())

ax.coastlines(linewidth=0.5)
plt.colorbar(p1)
ax.set_title("PI precipitation")

Text(0.5, 1.0, 'PI precipitation')
_images/1ac957e37fc3c50a7bf1a9a906130ba3e500e177b4076f5312eab5d1303377c3.png

Let’s use add_cyclic_point to get rid of the “white strip” in the plot#

fig, ax = plt.subplots(figsize=(3, 1.5), subplot_kw={
    'projection': ccrs.Robinson(central_longitude=210)})

# Note the differences in the next two lines
prec_PI_new, lon_new = add_cyclic_point(prec_PI, prec_PI.lon)
p1 = ax.contourf(lon_new, prec_PI.lat, prec_PI_new,
                 cmap='YlGnBu',
                 levels=np.linspace(0, 17, 18),
                 extend='both',
                 transform=ccrs.PlateCarree())

ax.coastlines(linewidth=0.5)
plt.colorbar(p1)
ax.set_title("PI precipitation")

Text(0.5, 1.0, 'PI precipitation')
_images/6ce464731e2ca828549f79105c6747faf4543d4b74a00d10730f75f14bf40b8a.png

Add your own plot of the midHolocene precipitation#

fig, ax = plt.subplots(figsize=(3, 1.5), subplot_kw={
    'projection': ccrs.Robinson(central_longitude=210)})
_images/706e90dd4d0df98925d4a0a3fcc2a6f105bf3bc9e645c9afba34a3d5189199a6.png

Add your own plot of the midHolocene - piControl#

fig, ax = plt.subplots(figsize=(3, 1.5), subplot_kw={
    'projection': ccrs.Robinson(central_longitude=210)})
_images/706e90dd4d0df98925d4a0a3fcc2a6f105bf3bc9e645c9afba34a3d5189199a6.png
Click here for the solution Copy and paste the code into the above cell 
prec_MH_new, lon_new = add_cyclic_point(prec_MH, prec_MH.lon)
p1 = ax.contourf(lon_new, prec_MH.lat, prec_MH_new,
                 cmap='YlGnBu',
                 levels=np.linspace(0, 17, 18),
                 extend='both',
                 transform=ccrs.PlateCarree())

ax.coastlines(linewidth=0.5)
plt.colorbar(p1)
ax.set_title("MH precipitation")


p1 = ax.contourf(lon_new, prec_MH.lat, prec_MH_new - prec_PI_new,
                 cmap='BrBG',
                 levels=np.linspace(-5, 5, 21),
                 extend='both',
                 transform=ccrs.PlateCarree())

ax.coastlines(linewidth=0.5)
plt.colorbar(p1)
ax.set_title("MH precipitation anomalies")

Discussion#

  • Do you see the ITCZ in the piControl?

  • How does the ITCZ change in the mid-Holocene?

  • Do you see changes of the monsoon precipitation? Is it consistent with findings from Kutzbach and Otto-Bliesner (1981, 1982)?

  • Ask the instructors about ITCZ and monsoon!

Analysis 4: how does the MH orbital forcing impact the atmospheric circulation?#

  • We plot the zonal mean zonal wind of PI in the Northern Hemisphere summer (JJA)

  • We use plt.gca().invert_yaxis() to invert the y-axis such that the high pressure is at the bottom

ds_PI.U

<xarray.DataArray 'U' (time: 12, lev: 26, lat: 96, lon: 144)> Size: 17MB
dask.array<concatenate, shape=(12, 26, 96, 144), dtype=float32, chunksize=(1, 26, 96, 144), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
  * lev      (lev) float64 208B 3.545 7.389 13.97 23.94 ... 929.6 970.6 992.6
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    mdims:         1
    units:         m/s
    long_name:     Zonal wind
    cell_methods:  time: mean

Plot the zonal mean of u-winds of JJA#

U_PI_za_jja = ds_PI.U.isel(time=slice(5, 8)).mean(('lon', 'time'))

U_PI_za_jja.plot.contourf(figsize=(3, 1.5),
                          levels=np.linspace(-50, 50, 21), extend='both', )
plt.gca().invert_yaxis()
_images/37ffd4e4b5c1c02dd5faf387eae101cd034a16c9bca454cba6db6cf0b9189169.png

Add your own plot to show the changes during the midHolocene#

Click here for the solution
Copy and paste the code into the above cell 
dU_MH_za_jja = (ds_MH.U - ds_PI.U).isel(time=slice(5, 8)).mean(('lon', 'time'))

dU_MH_za_jja.plot.contourf(figsize=(3, 1.5),
                          levels=np.linspace(-10, 10, 21), extend='both', )
plt.gca().invert_yaxis()

Did you notice any problem with the above plots?#

  • Check out the name of the vertical levels. It is called the hybrid sigma-pressure coordinate, which is NOT the pressure coordinate.

Atmosphere vertical coordinates

Figure: Atmosphere vertical coordinates

  • We need to interpolate from the hybrid sigma-pressure coordinate into the normal pressure coordinate

ds_PI.lev.standard_name

'atmosphere_hybrid_sigma_pressure_coordinate'

Use geocat interpolation to interpolate from hybrid sigma-pressure to pressure coordinate#

ds_PI.U

<xarray.DataArray 'U' (time: 12, lev: 26, lat: 96, lon: 144)> Size: 17MB
dask.array<concatenate, shape=(12, 26, 96, 144), dtype=float32, chunksize=(1, 26, 96, 144), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
  * lev      (lev) float64 208B 3.545 7.389 13.97 23.94 ... 929.6 970.6 992.6
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    mdims:         1
    units:         m/s
    long_name:     Zonal wind
    cell_methods:  time: mean
P_new_mb = np.array([
    1, 10, 20, 50., 100., 200., 300., 400., 500., 600., 700., 800., 900., 1000.
])

p0_mb = 1000.

Up_PI = interpolation.interp_hybrid_to_pressure(
    ds_PI.U,
    ds_PI.PS,
    ds_PI.hyam,
    ds_PI.hybm,
    p0=p0_mb*100.,
    new_levels=P_new_mb*100.,
    extrapolate=False)

Up_MH = interpolation.interp_hybrid_to_pressure(
    ds_MH.U,
    ds_MH.PS,
    ds_MH.hyam,
    ds_MH.hybm,
    p0=p0_mb*100.,
    new_levels=P_new_mb*100.,
    extrapolate=False)

Plot the correct results!#

fig, axes = plt.subplots(nrows=1, ncols=2,
                         figsize=(8, 2),
                         constrained_layout=True)

Up_PI.isel(time=slice(5, 8)).mean(('lon', 'time')).plot.contourf(
    ax=axes[0], levels=np.linspace(-50, 50, 21), extend='both')

(Up_MH - Up_PI).isel(time=slice(5, 8)).mean(('lon', 'time')).plot.contourf(
    ax=axes[1], levels=np.linspace(-10, 10, 21), extend='both')

for ax in axes:
    ax.invert_yaxis()
_images/77eb58b44da8e884a294475e7f9e917837a9400ade162cf6a6ff42e4bf204fd1.png

Small group discussion#

  • Which hemisphere has a stronger jet stream? Why? (remember that the plots are for JJA, the NH summer)

  • Does the mid-Holocene orbital forcing shift the jet stream?

  • Ask the instructors about the atmosphere circulation!

Analysis 5: how about sea-surface temperature?#

  • Ocean data is in ocn/hist

  • SST is the top level of TEMP

hist_dir = '/ocn/hist/'

files_PI_ocn = glob.glob(storage_dir + case_PI + hist_dir + '*.pop.h.0001*')
files_MH_ocn = glob.glob(storage_dir + case_MH + hist_dir + '*.pop.h.0001*')
print(*files_PI_ocn, sep='\n')
print(*files_MH_ocn, sep='\n')

ds_PI_ocn = xr.open_mfdataset(
    files_PI_ocn,
    combine="by_coords",
    data_vars="minimal",
    coords="minimal",
    compat="override",
)

ds_MH_ocn = xr.open_mfdataset(
    files_MH_ocn,
    combine="by_coords",
    data_vars="minimal",
    coords="minimal",
    compat="override",
)

# Again, we need this fix to get the correct time, i.e., month 1 to month 12
ds_PI_ocn['time'] = ds_PI_ocn.time.get_index('time') - timedelta(days=15)
ds_MH_ocn['time'] = ds_MH_ocn.time.get_index('time') - timedelta(days=15)

ds_PI_ocn

Calculate annual mean SST and make a simple plot using Xarray#

sst_PI = ds_PI_ocn.TEMP.isel(z_t=0).mean('time')

sst_MH = ds_MH_ocn.TEMP.isel(z_t=0).mean('time')

sst_PI.plot(size=1.5)

  • What do nlat and nlon mean?

  • Where is Greenland? It is the “north pole” in the model!
    POP ocean grid

Figure: POP ocean grid

Regridding is needed!#

  • Use the Regridder from xesmf

  • Use the util.grid_global to create a 1x1 regular grid

  • Use linear interpolation

  • Use IPython magic command to time the runtime

%%time

ds_PI_ocn['lat'] = ds_PI_ocn.TLAT
ds_PI_ocn['lon'] = ds_PI_ocn.TLONG

regridder = xesmf.Regridder(
    ds_in=ds_PI_ocn,
    ds_out=xesmf.util.grid_global(1, 1, cf=True, lon1=360),
    method='bilinear',
    periodic=True)

sst_PI_1x1 = regridder(sst_PI)

sst_MH_1x1 = regridder(sst_MH)

sst_PI_1x1

lat = sst_PI_1x1.lat
lon = sst_PI_1x1.lon

fig, axes = plt.subplots(nrows=1, ncols=3,
                         figsize=(10, 2),
                         subplot_kw={'projection': ccrs.Robinson(central_longitude=210)},
                         constrained_layout=True)

ax = axes[0]
ax.set_title("PI SST, annual mean")
sst_PI_new, lon_new = add_cyclic_point(sst_PI_1x1, lon)
p0 = ax.contourf(lon_new, lat, sst_PI_new,
                 levels=np.linspace(-2, 30, 17),
                 cmap='inferno', extend='both',
                 transform=ccrs.PlateCarree())
plt.colorbar(p0, ax=ax)

ax = axes[1]
ax.set_title("MH SST, annual mean")
sst_MH_new, lon_new = add_cyclic_point(sst_MH_1x1, lon)
p1 = ax.contourf(lon_new, lat, sst_MH_new,
                 levels=np.linspace(-2, 30, 17),
                 cmap='inferno', extend='both',
                 transform=ccrs.PlateCarree())
plt.colorbar(p1, ax=ax)

ax = axes[2]
ax.set_title("MH-PI SST")
p2 = ax.contourf(lon_new, lat, sst_MH_new - sst_PI_new,
                 cmap='coolwarm',
                 levels=np.linspace(-2, 2, 21),
                 extend='both',
                 transform=ccrs.PlateCarree())
plt.colorbar(p2, ax=ax)

for ax in axes:
    ax.set_global()
    ax.coastlines(linewidth=0.5)

# We could use savefig to save the plot as pdf
# plt.savefig('SST_xy.PI_vs_MH.pdf', format='pdf', bbox_inches="tight")

Small group discussion#

  • Colder SSTs over lots of the regions?

  • What’s the timescale of upper ocean response?

  • We need to consider internal variability and test significance. See 4b_opt2_analyze_long_midHolocene.

Summary so far#

  • Load CESM output from multiple files with Xarray

  • Plot variables as zonal means and in map view

  • Interpolate the atmospheric vertical coordinate from hybrid to pressure when you need pressure levels

  • Ocean models use an irregular grid with the poles displaced onto land, so regridding is needed

  • Mid-Holocene orbital forcing is seasonal and regional, and it affects atmospheric circulation, precipitation, and other climate states

Contents