Regional diagnostics - first year

Created by Ivan Lima on Mon Oct 18 2021 15:41:05 -0400

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

Last updated on Mon Oct 18 15:42:28 2021

from pencil_cases import *
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
# datadir = '/mnt/chromeos/removable/WD/Data/Postproc/1DPOP_CESM2/'
start_yr = 1851 # use "real" year to keep pandas datetime functionality

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

def plot_tseries(region, varname):
    ds1 = xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.005.ts_{}.nc'.format(region)),
                          decode_times=False, chunks='auto')
    ds2 = xr.open_dataset(os.path.join(datadir, 'cesm2.B1850.f09_g17.1dpop.006.ts_{}.nc'.format(region)),
                          decode_times=False, chunks='auto')
    ds1, ds2  = fixtime(ds1, start_yr), fixtime(ds2, start_yr)
    # select first year
    ds1 = ds1.sel(time=slice('1851-01-01','1851-12-31'))
    ds2 = ds2.sel(time=slice('1851-01-01','1851-12-31'))
    # convert daily to monthy averages
    df = pd.concat([ds1[varname].to_series(), ds2[varname].to_series()], axis=1, keys=['1dpop.005','1dpop.006'])
    # compute pi control run stats
    ds_ctrl = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.ts_{}.nc'.format(region)),
                              decode_times=False, chunks='auto')
    mean = ds_ctrl[varname].mean(dim='time').values.item()
    std  = ds_ctrl[varname].std(dim='time').values.item()
    df['pi control'] = mean
    df['min'] = mean - std
    df['max'] = mean + std
    # create plots
    p_tseries = df[['1dpop.005','1dpop.006','pi control']].hvplot(aspect=1.3, frame_height=350, grid=True, shared_axes=False, xlabel='time',
                                                                  ylabel='{}'.format(ds1[varname].units), title=varname)
    p_std = df[['min','max']].hvplot.area(y='min', y2='max', alpha=0.2, shared_axes=False, aspect=1.3, frame_height=350, grid=True,
                                          color='orange')
    return p_tseries * p_std

    # return df.hvplot(aspect=1.3, frame_height=350, grid=True, shared_axes=False, xlabel='time',
    #                  ylabel='{}'.format(ds1[varname].units), title=varname)

ds = xr.open_dataset(os.path.join(datadir, '{}.ts_near_global.nc'.format(case_list[0])), decode_times=False, chunks='auto')

cesm2.B1850.f09_g17.1dpop.005: 1D + G-terms + SSS restoring (tau = 30 days) + EBM off
cesm2.B1850.f09_g17.1dpop.006: 1D + G-terms + SSS restoring (tau = 365 days) + EBM off

The yellow shaded areas in the plots below represent the pi control mean plus/minus one standard deviation: $\bar{x} \pm \sigma$

Near global (60S - 60N)

plist_ng = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ng += plot_tseries('near_global',vname)

plist_ng.cols(2).opts(shared_axes=False)

North Atlantic

plist_na = hv.Layout()
for vname in ds.data_vars.keys():
    plist_na += plot_tseries('north_atl',vname)
    
plist_na.cols(2).opts(shared_axes=False)

Tropical Pacific

plist_tp = hv.Layout()
for vname in ds.data_vars.keys():
    plist_tp += plot_tseries('trop_pac',vname)
    
plist_tp.cols(2).opts(shared_axes=False)

Labrador Sea

plist_ls = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ls += plot_tseries('lab_sea',vname)
    
plist_ls.cols(2).opts(shared_axes=False)

Arctic

plist_ar = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ar += plot_tseries('arctic',vname)
    
plist_ar.cols(2).opts(shared_axes=False)

Southern Ocean

plist_so = hv.Layout()
for vname in ds.data_vars.keys():
    plist_so += plot_tseries('socean',vname)
    
plist_so.cols(2).opts(shared_axes=False)