Created by Ivan Lima on Mon Apr 12 2021 14:37:37 -0400
%matplotlib inline
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os, datetime, warnings
from tqdm import tnrange, notebook
from cesm_utils import fixtime
print('Last updated on {}'.format(datetime.datetime.now().ctime()))
Last updated on Mon Jul 12 14:31:01 2021
import seaborn as sns
sns.set_theme()
sns.set_context('paper')
warnings.filterwarnings('ignore')
plt.rcParams['figure.dpi'] = 100
from pencil_cases import *
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
start_yr = 1851 # use "real" year to keep pandas datetime functionality
caselist = ['cesm2.B1850.f09_g17.1dpop.004','cesm2.B1850.f09_g17.1dpop.005','cesm2.B1850.f09_g17.1dpop.006']
cases, cases_upper = {}, {}
droplist = ['TEMP','SALT','UVEL','VVEL','ULAT','ULONG',
'PT_INTERIOR','S_INTERIOR','UVEL_INTERIOR','VVEL_INTERIOR']
for run in sorted(caselist):
ds = xr.open_dataset(os.path.join(datadir, '{}.surf.nc'.format(run)), decode_times=False,
drop_variables=droplist, chunks='auto')
ds = fixtime(ds, start_yr)
ds['time'] = ds.time.to_index().to_datetimeindex()
cases[run] = ds.squeeze('z_t', drop=True)
ds2 = xr.open_dataset(os.path.join(datadir, '{}.upper_200m.nc'.format(run)), decode_times=False,
drop_variables=['UVEL','VVEL','PT_INTERIOR','S_INTERIOR',
'UVEL_INTERIOR','VVEL_INTERIOR','ULONG','ULAT'], chunks='auto')
ds2 = fixtime(ds2, start_yr)
ds2['time'] = ds2.time.to_index().to_datetimeindex()
cases_upper[run] = ds2
# case_info = {
# 'cesm2.B1850.f09_g17.1dpop.001': '1D restoring run',
# 'cesm2.B1850.f09_g17.1dpop.002': '1D monthly G-terms'
# }
def region_mean(ds, vname, regnum):
"""Compute time series of area weighted means for given region"""
data = ds[vname].where(ds.REGION_MASK==regnum)
weights = ds.TAREA.where(ds.REGION_MASK==regnum)
return ((data * weights).sum(dim=['nlat','nlon']) / weights.sum()).to_series()
# Labrador Sea REGION_MASK = 8
last_case = list(cases.keys())[-1]
labsea, units = {}, {}
for vname in notebook.tqdm(['HMXL','HBLT','IFRAC']):
datadict = {case_info[k]:region_mean(cases[k], vname, 8) for k in sorted(cases.keys())}
labsea[vname] = pd.DataFrame(datadict)
units[vname] = cases[last_case][vname].units
for vname in notebook.tqdm(['TEMP','SALT','RHO']):
datadict = {case_info[k]:region_mean(cases_upper[k], vname, 8) for k in sorted(cases_upper.keys())}
labsea[vname] = pd.DataFrame(datadict)
units[vname] = cases_upper[last_case][vname].units
# Arctic REGION_MASK = 10
arctic = {}
for vname in notebook.tqdm(['HMXL','HBLT','IFRAC']):
datadict = {case_info[k]:region_mean(cases[k], vname, 10) for k in sorted(cases.keys())}
arctic[vname] = pd.DataFrame(datadict)
for vname in notebook.tqdm(['TEMP','SALT','RHO']):
datadict = {case_info[k]:region_mean(cases_upper[k], vname, 10) for k in sorted(cases_upper.keys())}
arctic[vname] = pd.DataFrame(datadict)
def socean_mean(ds, vname):
"""Compute time series of area weighted means for Southern Ocean"""
data = ds[vname].where(ds.TLAT<=-60)
weights = ds.TAREA.where((ds.TLAT<=-60) & (ds.REGION_MASK!=0))
return ((data * weights).sum(dim=['nlat','nlon']) / weights.sum()).to_series()
# Southern Ocean
socean = {}
for vname in notebook.tqdm(['HMXL','HBLT','IFRAC']):
datadict = {case_info[k]:socean_mean(cases[k], vname) for k in sorted(cases.keys())}
socean[vname] = pd.DataFrame(datadict)
for vname in notebook.tqdm(['TEMP','SALT','RHO']):
datadict = {case_info[k]:socean_mean(cases_upper[k], vname) for k in sorted(cases_upper.keys())}
socean[vname] = pd.DataFrame(datadict)
The time series for TEMP, SALT and RHO represent the mean for the upper 200 m.
fig, axs = plt.subplots(3, 2, sharex=True, figsize=(11,11))
for ax, vname in zip(axs.ravel(), ['HMXL','HBLT','IFRAC','TEMP','SALT','RHO']):
_ = labsea[vname].plot(ax=ax, legend=False)
_ = ax.set(title='Labrador Sea', ylabel='{} ({})'.format(vname, units[vname]))
_ = ax.legend(loc='lower left',bbox_to_anchor=(1, 0))
fig, axs = plt.subplots(3, 2, sharex=True, figsize=(11,11))
for ax, vname in zip(axs.ravel(), ['HMXL','HBLT','IFRAC','TEMP','SALT','RHO']):
_ = arctic[vname].plot(ax=ax, legend=False)
_ = ax.set(title='Arctic', ylabel='{} ({})'.format(vname, units[vname]))
_ = ax.legend(loc='lower left',bbox_to_anchor=(1, 0))
fig, axs = plt.subplots(3, 2, sharex=True, figsize=(11,11))
for ax, vname in zip(axs.ravel(), ['HMXL','HBLT','IFRAC','TEMP','SALT','RHO']):
_ = socean[vname].plot(ax=ax, legend=False)
_ = ax.set(title='Southern Ocean', ylabel='{} ({})'.format(vname, units[vname]))
_ = ax.legend(loc='lower left',bbox_to_anchor=(1, 0))
Salinity mean in the first month
fig, axs = plt.subplots(1, 2, sharex=True, figsize=(9,4))
_ = socean['SALT'].loc[:'1851-02-01',:].plot(ax=axs[0], legend=False)
_ = socean['IFRAC'].loc[:'1851-02-01',:].plot(ax=axs[1], legend=False)
_ = axs[0].set(title='Southern Ocean', ylabel='{} ({})'.format('SALT', units['SALT']))
_ = axs[1].set(title='Southern Ocean', ylabel='{} ({})'.format('IFRAC', units['IFRAC']))
_ = axs[1].legend(loc='lower left',bbox_to_anchor=(1, 0))
The time series for TEMP, SALT and RHO represent the mean for the upper 200 m. The shaded area corresponds to one standard deviation from the climatological mean for the 3D run.
labsea_clim_3D = {}
varlist = ['HMXL','HBLT','IFRAC','TEMP','SALT','RHO']
for vname in varlist:
df = pd.read_csv('stats/labsea_{}.csv'.format(vname), index_col=0)
df.index.name = 'time'
df = df.rename(columns={'mean':'3D'})
labsea_clim_3D[vname] = df
labsea_clim = {}
for vname in labsea.keys():
labsea_clim[vname] = labsea[vname].groupby(labsea[vname].index.month).mean()
fig, axs = plt.subplots(3, 2, sharex=True, figsize=(11,11))
for ax, vname in zip(axs.ravel(), ['HMXL','HBLT','IFRAC','TEMP','SALT','RHO']):
_ = labsea_clim[vname].plot(ax=ax, legend=False)
_ = ax.set(title='Labrador Sea', xlabel='month', ylabel='{} ({})'.format(vname, units[vname]))
_ = labsea_clim_3D[vname]['3D'].plot(ax=ax, legend=legend)
_ = ax.fill_between(labsea_clim_3D[vname].index, labsea_clim_3D[vname]['3D']-labsea_clim_3D[vname]['std'],
labsea_clim_3D[vname]['3D']+labsea_clim_3D[vname]['std'], color='C3', alpha=0.3)
_ = ax.legend(loc='lower left',bbox_to_anchor=(1, 0))
arctic_clim_3D = {}
for vname in varlist:
df = pd.read_csv('stats/arctic_{}.csv'.format(vname), index_col=0)
df.index.name = 'time'
df = df.rename(columns={'mean':'3D'})
arctic_clim_3D[vname] = df
arctic_clim = {}
for vname in arctic.keys():
arctic_clim[vname] = arctic[vname].groupby(arctic[vname].index.month).mean()
fig, axs = plt.subplots(3, 2, sharex=True, figsize=(11,11))
for ax, vname in zip(axs.ravel(), ['HMXL','HBLT','IFRAC','TEMP','SALT','RHO']):
_ = arctic_clim[vname].plot(ax=ax, legend=False)
_ = ax.set(title='Arctic', xlabel='month', ylabel='{} ({})'.format(vname, units[vname]))
_ = arctic_clim_3D[vname]['3D'].plot(ax=ax, legend=legend)
_ = ax.fill_between(arctic_clim_3D[vname].index, arctic_clim_3D[vname]['3D']-arctic_clim_3D[vname]['std'],
arctic_clim_3D[vname]['3D']+arctic_clim_3D[vname]['std'], color='C3', alpha=0.3)
_ = ax.legend(loc='lower left',bbox_to_anchor=(1, 0))
socean_clim_3D = {}
for vname in varlist:
df = pd.read_csv('stats/socean_{}.csv'.format(vname), index_col=0)
df.index.name = 'time'
df = df.rename(columns={'mean':'3D'})
socean_clim_3D[vname] = df
socean_clim = {}
for vname in socean.keys():
socean_clim[vname] = socean[vname].groupby(socean[vname].index.month).mean()
fig, axs = plt.subplots(3, 2, sharex=True, figsize=(11,11))
for ax, vname in zip(axs.ravel(), ['HMXL','HBLT','IFRAC','TEMP','SALT','RHO']):
_ = socean_clim[vname].plot(ax=ax, legend=False)
_ = ax.set(title='Southern Ocean', xlabel='month', ylabel='{} ({})'.format(vname, units[vname]))
_ = socean_clim_3D[vname]['3D'].plot(ax=ax, legend=legend)
_ = ax.fill_between(socean_clim_3D[vname].index, socean_clim_3D[vname]['3D']-socean_clim_3D[vname]['std'],
socean_clim_3D[vname]['3D']+socean_clim_3D[vname]['std'], color='C3', alpha=0.3)
_ = ax.legend(loc='lower left',bbox_to_anchor=(1, 0))
