Restoring runs¶

Created by Ivan Lima on Tue Jul 12 2022 13:31:50 -0400

In [1]:
import xarray as xr
import pandas as pd
import numpy as np
import holoviews as hv
import hvplot.xarray
import hvplot.pandas
import os, datetime, warnings
import my_bokeh_themes
from cesm_utils import fixtime, spd
warnings.filterwarnings('ignore')
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Tue Jul 12 14:33:01 2022
In [2]:
from pencil_cases import *
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
start_yr = 1851

case_list = ['cesm2.B1850.f09_g17.1dpop.010', 'cesm2.B1850.f09_g17.1dpop.015']
for case in sorted(case_list):
    print('{}: {}'.format(case, case_info[case]))

def get_weights(ds, vname):
    """Get weights & dimensions to compute weighted means"""
    if 'z_t' in ds[vname].dims:
        weights = ds.TAREA * ds.dz
        dims = ['nlat','nlon','z_t']
    else:
        weights = ds.TAREA
        dims = ['nlat','nlon']

    weights = weights.where(ds[vname].isel(time=0).notnull())
    return weights, dims

def global_mean(ds, vname, upper200m=False, surf=False):
    if upper200m:
        ds = ds.isel(z_t=slice(0,20)) # upper 200 m
    elif surf:
        ds = ds.isel(z_t=0) # surface
    weights, dims = get_weights(ds, vname)
    da = ds[vname].weighted(weights.fillna(0)).mean(dim=dims)
    da.attrs = ds[vname].attrs
    da.time.attrs = ds.time.attrs
    return da

def plot_tseries(vname, upper200m=False, surf=False):
    data_dict = {}
    cols = []
    for case in case_list:
        col_name = '.'.join(case.split('.')[-2:])
        cols.append(col_name)
        with xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case, vname)), decode_times=False, chunks={'time':10}) as ds:
            ds = fixtime(ds, start_yr)
            ds['time'] = ds.time.to_index().to_datetimeindex().year
            units = ds[vname].units
            data_dict[col_name] = global_mean(ds, vname, upper200m, surf).to_series()

    df = pd.DataFrame(data_dict)
    df.columns.name = 'Case'
    if upper200m:
        title='Global mean {} (upper 200 m)'.format(vname)
        ds_ctrl = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.ts_upper200m.nc'),
                                  decode_times=False, chunks='auto')
    elif surf:
        title='Global mean {} (surface)'.format(vname)
        ds_ctrl = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.ts_surf.nc'),
                                  decode_times=False, chunks='auto')
    else:
        title='Global mean {}'.format(vname)
        ds_ctrl = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.ts_global.nc'),
                                  decode_times=False, chunks='auto')

    mean = ds_ctrl[vname].mean(dim='time').values.item()
    std  = ds_ctrl[vname].std(dim='time').values.item()
    df['lo'] = mean - std
    df['hi'] = mean + std
    df['min'] = ds_ctrl[vname].min(dim='time').values.item()
    df['max'] = ds_ctrl[vname].max(dim='time').values.item()
    
    p_tseries = df[cols].hvplot(aspect=1.1, frame_height=350, grid=True, title=title, ylabel=units)
    p_std = df[['lo','hi']].hvplot.area(y='lo', y2='hi', alpha=0.2, frame_height=350, grid=True, color='green')
    p_minmax = df[['min','max']].hvplot.area(y='min', y2='max', alpha=0.2, frame_height=350, grid=True, color='lightgreen')
    
    return p_tseries * p_std * p_minmax

cesm2.B1850.f09_g17.1dpop.010: 1D + T, S, U & V restoring (tau = 15 days) + EBM off + start year = 1161
cesm2.B1850.f09_g17.1dpop.015: 1D + T, S, U & V restoring (tau = 30 days) + EBM off + start year = 1161

Global mean time series¶

The green shaded areas in the plots below represent the pi control mean plus/minus one standard deviation: $\bar{x} \pm \sigma$. The light green shaded area represents the pi control min and max range.

In [3]:
plist = hv.Layout()
for vname in ['TEMP','SALT']:
    plist += plot_tseries(vname)

plist.cols(2).opts(shared_axes=False)
Out[3]:

Global mean time series (upper 200 m)¶

In [4]:
plist = hv.Layout()
for vname in ['TEMP','SALT']:
    plist += plot_tseries(vname, upper200m=True)

plist.cols(2).opts(shared_axes=False)
Out[4]:

Global mean time series (surface)¶

In [5]:
plist = hv.Layout()
for vname in ['TEMP','SALT']:
    plist += plot_tseries(vname, surf=True)

plist.cols(2).opts(shared_axes=False)
Out[5]:

Global mean profile time series¶

In [6]:
rest_units = {'PT_INTERIOR': 'deg C/century', 'S_INTERIOR': 'msu/century'}

def plot_profile_tseries(case, vname, anom=True):
    with xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case, vname)), 
                         decode_times=False, chunks={'time':10}) as ds:
        ds = fixtime(ds, start_yr)
        ds['time'] = ds.time.to_index().to_datetimeindex().year
        w = ds.TAREA.where(ds[vname].isel(time=0, z_t=0).notnull())
        if vname in ['PT_INTERIOR', 'S_INTERIOR']:
            da = ds[vname].weighted(w.fillna(0)).mean(dim=['nlat','nlon'])  * spd * 365 * 100 # X/sec -> X/century
            units = rest_units[vname]
        else:
            da = ds[vname].weighted(w.fillna(0)).mean(dim=['nlat','nlon'])
            units = ds[vname].units
        da['z_t'] = da.z_t/100 # cm -> m
        case_name = '.'.join(case.split('.')[-2:])        
        if anom:
            da = da - da.mean(dim='time') # anomaly from mean global profile
            vmax = np.abs([da.values.min(), da.values.max()]).max()
            clim=(-vmax, vmax)
        else:
            # clim=(None,None)
            # cmap='kbc'
            vmax = np.abs([da.values.min(), da.values.max()]).max()
            clim=(-vmax, vmax)
        return da.hvplot.quadmesh(x='time', y='z_t', z=vname, cmap='bwr', frame_height=350, aspect=1.2, clim=clim,
                                  title='Global {} profile ({})'.format(vname, case_name), xaxis='top',
                                  clabel='{} ({})'.format(vname, units)).opts(invert_yaxis=True)

def plot_mean_profile(case, vname):
    with xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case, vname)), 
                         decode_times=False, chunks={'time':10}) as ds:
        w = ds.TAREA.where(ds[vname].isel(time=0, z_t=0).notnull())
        da = ds[vname].weighted(w.fillna(0)).mean(dim=['time','nlat','nlon'])  * spd * 365 * 100 # X/sec -> X/century
        da['z_t'] = da.z_t/100 # cm -> m
        case_name = '.'.join(case.split('.')[-2:])
        return da.hvplot.line(invert=True, frame_height=350, frame_width=350, ylabel=units[vname], xaxis='top',
                              title='{} ({})'.format(vname,case_name), grid=True, label='full time series').opts(invert_yaxis=True)

The plots below show anomalies from the mean global profile.

Temperature¶

In [7]:
plist_prof = hv.Layout()
for case in case_list:
    plist_prof += plot_profile_tseries(case, 'TEMP')

plist_prof.cols(2).opts(shared_axes=False)
Out[7]:

Salinity¶

In [8]:
plist_prof = hv.Layout()
for case in case_list:
    plist_prof += plot_profile_tseries(case, 'SALT')

plist_prof.cols(2).opts(shared_axes=False)
Out[8]:

Restoring fluxes (g-terms) for TEMP & SALT from "restoring" runs¶

In [9]:
ptemp = plot_profile_tseries('cesm2.B1850.f09_g17.1dpop.010', 'PT_INTERIOR', anom=False)
psalt = plot_profile_tseries('cesm2.B1850.f09_g17.1dpop.010', 'S_INTERIOR', anom=False)
(ptemp + psalt).cols(2).opts(shared_axes=False)
Out[9]:
In [10]:
ptemp = plot_profile_tseries('cesm2.B1850.f09_g17.1dpop.015', 'PT_INTERIOR', anom=False)
psalt = plot_profile_tseries('cesm2.B1850.f09_g17.1dpop.015', 'S_INTERIOR', anom=False)
(ptemp + psalt).cols(2).opts(shared_axes=False)
Out[10]:

Mean restoring fluxes (g-terms) profiles for TEMP & SALT¶

In [11]:
units = {'PT_INTERIOR': 'deg C/century', 'S_INTERIOR': 'msu/century'}
pt = plot_mean_profile('cesm2.B1850.f09_g17.1dpop.010', 'PT_INTERIOR')
ps = plot_mean_profile('cesm2.B1850.f09_g17.1dpop.010', 'S_INTERIOR')
(pt + ps).cols(2).opts(shared_axes=False)
Out[11]:
In [12]:
pt = plot_mean_profile('cesm2.B1850.f09_g17.1dpop.015', 'PT_INTERIOR')
ps = plot_mean_profile('cesm2.B1850.f09_g17.1dpop.015', 'S_INTERIOR')
(pt + ps).cols(2).opts(shared_axes=False)
Out[12]:

Mean restoring fluxes (g-terms) for TEMP & SALT¶

In [13]:
data_dict, means = {}, {}
for vname in ['PT_INTERIOR', 'S_INTERIOR']:
    with xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.010.{}.annual.nc'.format(vname)),
                         decode_times=False, chunks={'time':10}) as ds:
        ds = fixtime(ds, start_yr)
        ds['time'] = ds.time.to_index().to_datetimeindex().year
        data_dict[vname] = global_mean(ds, vname).to_series() * spd * 365 * 100 # X/sec -> X/century
        means[vname] = data_dict[vname].mean()

df = pd.DataFrame(data_dict)
plist = df.hvplot(aspect=1.1, frame_height=350, grid=True, subplots=True, shared_axes=False)
for p in plist.items():
    p[1].opts(title='Global mean {} (mean={:.4g})'.format(p[0], means[p[0]]), ylabel=units[p[0]])
plist
Out[13]:
In [14]:
data_dict, means = {}, {}
for vname in ['PT_INTERIOR', 'S_INTERIOR']:
    with xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.015.{}.annual.nc'.format(vname)),
                         decode_times=False, chunks={'time':10}) as ds:
        ds = fixtime(ds, start_yr)
        ds['time'] = ds.time.to_index().to_datetimeindex().year
        data_dict[vname] = global_mean(ds, vname).to_series() * spd * 365 * 100 # X/sec -> X/century
        means[vname] = data_dict[vname].mean()

df = pd.DataFrame(data_dict)
plist = df.hvplot(aspect=1.1, frame_height=350, grid=True, subplots=True, shared_axes=False)
for p in plist.items():
    p[1].opts(title='Global mean {} (mean={:.4g})'.format(p[0], means[p[0]]), ylabel=units[p[0]])
plist
Out[14]: