Compare initial fields with restart fields

Created by Ivan Lima on Tue Apr 27 2021 09:10:41 -0400

%matplotlib inline
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from cartopy import crs as ccrs
from cartopy import feature as cfeature
import os, datetime, warnings, cmocean
from matplotlib import colors
from cesm_utils import fixtime, spd, da2ma
from mpl_utils import center_cmap
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Tue Apr 27 11:21:31 2021

warnings.filterwarnings('ignore')
plt.rcParams['figure.dpi'] = 100

Model output

datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2/init/'

ds_initial = xr.open_dataset(datadir + 'cesm2.B1850.f09_g17.1dpop.002.pop.h.0001-01-01.nc',
                             decode_times=False)
ds_initial = ds_initial.squeeze('time', drop=True)

salt_model = ds_initial.SALT[:20,:,:].weighted(ds_initial.dz[:20]).mean('z_t')
salt_model = salt_model.where(ds_initial.REGION_MASK>0)
tlat, tlon = ds_initial.TLAT, ds_initial.TLONG
tlat, tlon = np.c_[tlat,tlat[:,0]], np.c_[tlon,tlon[:,0]]

fig, ax = plt.subplots(figsize=(10,10), subplot_kw={'projection':ccrs.Robinson(central_longitude=210)})
_ = ax.set_global()
_ = ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
_ = ax.coastlines(linewidth=0.5)
_ = ax.set_title('model')
pm = ax.pcolormesh(tlon, tlat, salt_model, transform=ccrs.PlateCarree(), cmap=cmocean.cm.haline)
cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.025)
cb.set_label('average salinity in upper 200 m (g/kg)')

Restart file

ds_restart = xr.open_dataset(datadir + 'b.e21.B1850.f09_g17.CMIP6-piControl.001.pop.r.0601-01-01-00000.nc',
                             decode_times=False)
data = ds_restart.SALT_CUR.values * 1000
marr = da2ma(ds_initial.SALT)
data = np.ma.masked_where(marr.mask, data)
weights = ds_initial.dz.values

salt_restart = np.ma.average(data[:20,:,:], weights=weights[:20], axis=0)
salt_restart = np.ma.masked_where(ds_initial.REGION_MASK<0, salt_restart)

fig, ax = plt.subplots(figsize=(10,10), subplot_kw={'projection':ccrs.Robinson(central_longitude=210)})
_ = ax.set_global()
_ = ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
_ = ax.coastlines(linewidth=0.5)
_ = ax.set_title('restart file')
pm = ax.pcolormesh(tlon, tlat, salt_restart, transform=ccrs.PlateCarree(), cmap=cmocean.cm.haline)
cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.025)
cb.set_label('average salinity in upper 200 m (g/kg)')

Difference

diff = da2ma(salt_model) - salt_restart

fig, ax = plt.subplots(figsize=(10,10), subplot_kw={'projection':ccrs.Robinson(central_longitude=210)})
_ = ax.set_global()
_ = ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
_ = ax.coastlines(linewidth=0.5)
_ = ax.set_title('restart file')
pm = ax.pcolormesh(tlon, tlat, diff, transform=ccrs.PlateCarree(), cmap=cmocean.cm.balance, norm=colors.TwoSlopeNorm(0))
# center_cmap(pm)
cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.025)
cb.set_label('average salinity in upper 200 m (g/kg)')