3. Analyze CESM Output

3. Analyze CESM Output#

Tutorials at the 2025 paleoCAMP | June 16–June 30, 2025

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 temperature, precipitation, and sea-surface temperature

Time to learn: 60 minutes

How to get started?

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

Start the JupyterHub. From the dropdown menu, select

Casper PBS Batch for Resource Selection
casper for Queue or Reservation
UAZN0042 for Project Account
06:00:00 for Wall Time

Jupyter Login

Figure: Jupyterhub Casper

Lauch a terminal: click File, New, and Terminal

Launch Terminal

Figure: Launch a Terminal

Type git clone https://github.com/jiang-zhu/paleocamp2025.git
Find and open 3_analyze_CESM_output.ipynb from the left sidebar
Make sure to use NPL 2025a as your Kernel (select from top right corner)

Figure: Launch a Terminal

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

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_quick_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 month 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)
ds_MH = xr.open_mfdataset(files_MH)

ds_PI

<xarray.Dataset> Size: 541MB
Dimensions:       (time: 12, zlon: 1, nbnd: 2, lat: 96, lev: 26, ilev: 27,
                   lon: 144)
Coordinates:
  * lat           (lat) float64 768B -90.0 -88.11 -86.21 ... 86.21 88.11 90.0
  * zlon          (zlon) float64 8B 0.0
  * lon           (lon) float64 1kB 0.0 2.5 5.0 7.5 ... 350.0 352.5 355.0 357.5
  * 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
Dimensions without coordinates: nbnd
Data variables: (12/122)
    zlon_bnds     (time, zlon, nbnd) float64 192B dask.array<chunksize=(1, 1, 2), meta=np.ndarray>
    gw            (time, lat) float64 9kB dask.array<chunksize=(1, 96), meta=np.ndarray>
    hyam          (time, lev) float64 2kB dask.array<chunksize=(1, 26), meta=np.ndarray>
    hybm          (time, lev) float64 2kB dask.array<chunksize=(1, 26), meta=np.ndarray>
    P0            (time) float64 96B 1e+05 1e+05 1e+05 ... 1e+05 1e+05 1e+05
    hyai          (time, ilev) float64 3kB dask.array<chunksize=(1, 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:              derecho2
    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

# We need this fix 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:       (time: 12, zlon: 1, nbnd: 2, lat: 96, lev: 26, ilev: 27,
                   lon: 144)
Coordinates:
  * lat           (lat) float64 768B -90.0 -88.11 -86.21 ... 86.21 88.11 90.0
  * zlon          (zlon) float64 8B 0.0
  * lon           (lon) float64 1kB 0.0 2.5 5.0 7.5 ... 350.0 352.5 355.0 357.5
  * 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
Dimensions without coordinates: nbnd
Data variables: (12/122)
    zlon_bnds     (time, zlon, nbnd) float64 192B dask.array<chunksize=(1, 1, 2), meta=np.ndarray>
    gw            (time, lat) float64 9kB dask.array<chunksize=(1, 96), meta=np.ndarray>
    hyam          (time, lev) float64 2kB dask.array<chunksize=(1, 26), meta=np.ndarray>
    hybm          (time, lev) float64 2kB dask.array<chunksize=(1, 26), meta=np.ndarray>
    P0            (time) float64 96B 1e+05 1e+05 1e+05 ... 1e+05 1e+05 1e+05
    hyai          (time, ilev) float64 3kB dask.array<chunksize=(1, 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:              derecho2
    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 signle 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:
  * 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
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
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)

/glade/u/apps/opt/conda/envs/npl-2025a/lib/python3.12/site-packages/dask/config.py:787: FutureWarning: Dask configuration key 'allowed-failures' has been deprecated; please use 'distributed.scheduler.allowed-failures' instead
  warnings.warn(

<matplotlib.collections.QuadMesh at 0x1489712bfbc0>

_images/48359b936f0cdc7a23813385b97e8624f0c3bc8e2d1868c838abd4211f49d238.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:
  * 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
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
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/2532300d196468d65cccac37253fe73c12b819b7a69dc3fdcc0017a90bb57e38.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 midHolocene?#

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 degree

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

xarray.DataArray

'lat'

lat: 96

-90.0 -88.11 -86.21 -84.32 -82.42 ... 82.42 84.32 86.21 88.11 90.0

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: (1)

lat

(lat)

float64

-90.0 -88.11 -86.21 ... 88.11 90.0

long_name :: latitude
units :: degrees_north

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.      ])

Indexes: (1)

lat

PandasIndex

PandasIndex(Index([              -90.0,  -88.10526315789474,  -86.21052631578948,
         -84.3157894736842,  -82.42105263157895,  -80.52631578947368,
        -78.63157894736842,  -76.73684210526316,  -74.84210526315789,
        -72.94736842105263,  -71.05263157894737,  -69.15789473684211,
        -67.26315789473685,  -65.36842105263158, -63.473684210526315,
       -61.578947368421055, -59.684210526315795,  -57.78947368421053,
        -55.89473684210527,               -54.0,  -52.10526315789474,
        -50.21052631578947,  -48.31578947368421,  -46.42105263157895,
       -44.526315789473685, -42.631578947368425,  -40.73684210526316,
         -38.8421052631579,  -36.94736842105264,  -35.05263157894737,
        -33.15789473684211, -31.263157894736842, -29.368421052631582,
       -27.473684210526322, -25.578947368421055, -23.684210526315795,
       -21.789473684210535,  -19.89473684210526,               -18.0,
        -16.10526315789474,  -14.21052631578948,  -12.31578947368422,
       -10.421052631578945,  -8.526315789473685,  -6.631578947368425,
        -4.736842105263165, -2.8421052631579045, -0.9473684210526301,
        0.9473684210526301,  2.8421052631578902,    4.73684210526315,
         6.631578947368411,   8.526315789473685,  10.421052631578945,
        12.315789473684205,  14.210526315789465,  16.105263157894726,
                      18.0,   19.89473684210526,   21.78947368421052,
         23.68421052631578,   25.57894736842104,  27.473684210526315,
        29.368421052631575,  31.263157894736835,  33.157894736842096,
        35.052631578947356,   36.94736842105263,   38.84210526315789,
         40.73684210526315,   42.63157894736841,   44.52631578947367,
         46.42105263157893,   48.31578947368419,   50.21052631578948,
         52.10526315789474,                54.0,   55.89473684210526,
         57.78947368421052,   59.68421052631578,   61.57894736842104,
          63.4736842105263,   65.36842105263156,   67.26315789473682,
         69.15789473684211,   71.05263157894737,   72.94736842105263,
         74.84210526315789,   76.73684210526315,   78.63157894736841,
         80.52631578947367,   82.42105263157893,   84.31578947368419,
         86.21052631578945,   88.10526315789474,                90.0],
      dtype='float64', name='lat'))

Attributes: (2)
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.3271967819253 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 mid-Holocene#

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#

Should the mid-Holocene be warmer or colder than the preindustrial?
- Geological records suggest that MH may be warmer than the present day
- Climate models in general suggest a colder mid-Holocene
- Further reading: Osman et al. (2021)
- Ask Jess and Jiang about the Holocene Temperature Conundrum

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 Agugust 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:
  * 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
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
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:
  * 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
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
Attributes:
    units:         m/s
    long_name:     Large-scale (stable) precipitation rate (liq + ice)
    cell_methods:  time: mean

m_p_s_to_mm_p_day = 86400000

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: 111kB
dask.array<mul, shape=(96, 144), dtype=float64, 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

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/e684b8536bcd1357f4bf3c9c30537cebf781fd7c9b1d61b0f6eaff63dad4a16d.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/2f852a091d7863c63b4cc9ffc8bd0aa02391118318a1cacd714d25f61185bdf1.png

Add your own plot of the midHolocene precipitatin#

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

_images/d631faa0a1a61ece7a84f826274ca56b5e8d7f2fc6e80d93b9ccd87bb4c6a59c.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)})

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 Kathleen, Jane, Daniel, and Kevin 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:
  * 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
  * lev      (lev) float64 208B 3.545 7.389 13.97 23.94 ... 929.6 970.6 992.6
  * time     (time) object 96B 0001-01-17 00:00:00 ... 0001-12-17 00:00:00
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/2a57fb1b77723f9405909d79934cb0a9f7f2570bd71f3e9f28567f615627ec65.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 hybrid sigma-pressure coordinate, which is NOT the pressure coordinate.

Atmopshere vertical coordinates

Figure: Atmopshere 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 presure coordinate#

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/8e91326b81ee51c94de0a67eac8f95bdae6b10830d6687782905848b68ab8fd5.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 Sylvia, Jane and Kevin 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)
ds_MH_ocn = xr.open_mfdataset(files_MH_ocn)

# 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)

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

ds_PI_ocn

<xarray.Dataset> Size: 20GB
Dimensions:                 (time: 12, d2: 2, z_t: 60, z_w: 60, nlat: 384,
                             nlon: 320, z_w_bot: 60, z_w_top: 60,
                             transport_reg: 2, moc_comp: 3, moc_z: 61,
                             lat_aux_grid: 395, transport_comp: 5, z_t_150m: 15)
Coordinates: (12/15)
    moc_components          (moc_comp) |S384 1kB dask.array<chunksize=(3,), meta=np.ndarray>
    transport_components    (transport_comp) |S384 2kB dask.array<chunksize=(5,), meta=np.ndarray>
    transport_regions       (transport_reg) |S384 768B dask.array<chunksize=(2,), meta=np.ndarray>
  * time                    (time) object 96B 0001-01-17 00:00:00 ... 0001-12...
  * z_t                     (z_t) float32 240B 500.0 1.5e+03 ... 5.375e+05
  * z_t_150m                (z_t_150m) float32 60B 500.0 1.5e+03 ... 1.45e+04
    ...                      ...
  * lat_aux_grid            (lat_aux_grid) float32 2kB -79.49 -78.95 ... 90.0
  * moc_z                   (moc_z) float32 244B 0.0 1e+03 ... 5.25e+05 5.5e+05
    ULONG                   (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
    ULAT                    (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
    TLONG                   (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
    TLAT                    (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
Dimensions without coordinates: d2, nlat, nlon, transport_reg, moc_comp,
                                transport_comp
Data variables: (12/171)
    time_bound              (time, d2) object 192B dask.array<chunksize=(1, 2), meta=np.ndarray>
    dz                      (time, z_t) float32 3kB dask.array<chunksize=(1, 60), meta=np.ndarray>
    dzw                     (time, z_w) float32 3kB dask.array<chunksize=(1, 60), meta=np.ndarray>
    KMT                     (time, nlat, nlon) float64 12MB dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    KMU                     (time, nlat, nlon) float64 12MB dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    REGION_MASK             (time, nlat, nlon) float64 12MB dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    ...                      ...
    XBLT                    (time, nlat, nlon) float32 6MB dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    TBLT                    (time, nlat, nlon) float32 6MB dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    BSF                     (time, nlat, nlon) float32 6MB dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    MOC                     (time, transport_reg, moc_comp, moc_z, lat_aux_grid) float32 7MB dask.array<chunksize=(1, 2, 3, 61, 395), meta=np.ndarray>
    N_HEAT                  (time, transport_reg, transport_comp, lat_aux_grid) float32 190kB dask.array<chunksize=(1, 2, 5, 395), meta=np.ndarray>
    N_SALT                  (time, transport_reg, transport_comp, lat_aux_grid) float32 190kB dask.array<chunksize=(1, 2, 5, 395), meta=np.ndarray>
Attributes:
    title:             b.e21.B1850.f19_g17.piControl.001
    history:           none
    Conventions:       CF-1.0; http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-cu...
    time_period_freq:  month_1
    model_doi_url:     https://doi.org/10.5065/D67H1H0V
    contents:          Diagnostic and Prognostic Variables
    source:            CCSM POP2, the CCSM Ocean Component
    revision:          $Id$
    calendar:          All years have exactly  365 days.
    start_time:        This dataset was created on 2025-06-14 at 12:35:19.2
    cell_methods:      cell_methods = time: mean ==> the variable values are ...

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)

<matplotlib.collections.QuadMesh at 0x14895d798590>

_images/fa0755451b4895fb2ac9c75faf35d90c98dbbfa87ab74a8f57fff1932ee193bf.png

What do nlat and nlon mean?
Where is Greenland? It is the “north pole” in the model!

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

CPU times: user 4.79 s, sys: 279 ms, total: 5.07 s
Wall time: 5.64 s

<xarray.DataArray (lat: 180, lon: 360)> Size: 259kB
dask.array<astype, shape=(180, 360), dtype=float32, chunksize=(180, 360), chunktype=numpy.ndarray>
Coordinates:
    z_t                 float32 4B 500.0
  * lon                 (lon) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    latitude_longitude  float64 8B nan
  * lat                 (lat) float64 1kB -89.5 -88.5 -87.5 ... 87.5 88.5 89.5
Attributes:
    regrid_method:  bilinear

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")

_images/958ffa06e1ca9d9fdf87ea357bf89ba7e29da4c93623a26a3ed68bf82af3d37e.png

Small group discussion#

Colder SSTs over lots of the regions? Internal variability?
We need to consider internal variability and test significance
What’s the timescale of upper ocean response?

Summary so far#

Load multiple CESM output using Xarray
Plot variables in zonal mean and map view
Interpolate the atmosphere vertical coordinate from hybrid to pressure is needed!
Ocean models usually use irregular grid with poles moved into land; regridding is needed!
Mid-Holocene orbital forcing is seasonal and impacts atmosphere circulation, precipitation, and other climate states.

3. Analyze CESM Output

Contents

3. Analyze CESM Output#

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

Load data#

Compute the zonal mean of solar insolation and make plots#

Small group discussion#

Analysis 2: Does the Earth receive more radiation during the midHolocene?#

Example calculation for the piControl#

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

Add your own calculation for the mid-Holocene#

Small group discussion#

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

Plot total precipitation of piControl#

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

Add your own plot of the midHolocene precipitatin#

Add your own plot of the midHolocene - piControl#

Discussion#

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

Plot the zonal mean of u-winds of JJA#

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

Did you notice any problem with the above plots?#

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

Plot the correct results!#

Small group discussion#

Analysis 5: how about sea-surface temperature?#

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

Regridding is needed!#

Small group discussion#

Summary so far#

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