Surface tritium from model & observations

Created by Ivan Lima on Mon Apr 30 2018 14:13:07 -0400

In [1]:
%matplotlib notebook
import pandas as pd
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import os, gsw
from matplotlib import colors
from cesm_utils import da2ma, fixtime, spd, dpm
from he3_inventory import *
from cases import *
from datetime import datetime
print('Last updated on %s'%datetime.now().ctime())
plt.rcParams['xtick.minor.bottom'] = False # fix strange minor ticks in time series plots

Last updated on Tue Apr 30 09:43:03 2019

Case description

In [2]:
case = 'tr3he.GIAF.06'
print(case,':',case_info[case])

tr3he.GIAF.06 : 186-year run with interannualy varying forcing (CORE 1948-2009). Physics is
spun-up for 126 years (2 IAF cycles: 1824-1949). Model is spun-up with pre-bomb tritium fluxes and
historic tritium deposition starts in the third IAF cycle (1950-2009). Time & latitude dependent
tritium concentration in precipitation (Cp) is estimated from GNIP data (Michaela's blended fits).
Tritium concentration in water vapor (Cv) is set to 0.7 Cp. Relative humidity has a max value of
0.95. Added tritium river fluxes. Time & latitude dependent tritium concentration in rivers is
estimated from GNIP data (Michaela's blended fits). Tritium concentration in precipitation over
ocean and land are smoothed over a 1 degree latitudinal grid. Run includes vapor flux correction.
In [3]:
def get_lat_band(df, latmin, latmax, zsurf=-10):
    """ extract obs at given latitude bands """
    band = df[(df.tritium_flag==2) | (df.tritium_flag==6)] # only valid measurements
    band = band[band.latitude <= latmax]
    band = band[band.latitude > latmin]
    if zsurf:
        band = band[band.depth > zsurf]
    band = band[['stadate','latitude','longitude','depth','tritium','tritium_error','tritium_flag']]
    return band.sort_values('stadate') # sort by date

def lat_band_mean(ds, latmin, latmax, vname='TRITIUM'): 
    """ compute min, max & area weighted averages for model """
    da = ds[vname].where((ds.TLAT>latmin) & (ds.TLAT<=latmax))
    area = ds.TAREA.where((ds.TLAT>latmin) & (ds.TLAT<=latmax))
    vmean = (da * area).sum(dim=['nlat','nlon']) / area.sum(dim=['nlat','nlon'])
    vstd = np.sqrt(((da - vmean)**2 * area).sum(dim=['nlat','nlon']) / area.sum(dim=['nlat','nlon']))
    vmin, vmax = vmean - 1.96 * vstd, vmean + 1.96 * vstd # 95% confidence interval
    vmin[vmin<0] = 0
    return pd.DataFrame({'model':vmean.to_series(), 'vmin': vmin.to_series(), 'vmax':vmax.to_series()}) 

lat_bands = ((-90,-60), (-60,-45), (-45,-5), (-5,15), (15,45), (45,60), (60,90))
zbins2 = np.array([-4000, -3500, -3000, -2500, -2000, -1500, -1000, -500, -400, -300, -250, -200, -150, -100, -50, 0])

Read tritium database

In [4]:
def dyear2date(strn_list, yr0=1950): # convert decimal year to date objects
    ndays = [np.round((np.float(x)-yr0) * 365.2425) for x in strn_list] # number of days since 1950-01-01
    return pd.to_datetime(ndays, unit='D',origin='%d-01-01'%yr0)

database = pd.read_csv('../data/tritium-helium_database/tritium_helium_all_revised.csv', index_col=0, na_values=-999,
                       parse_dates=[6], date_parser=dyear2date, low_memory=False)

database['depth'] = gsw.z_from_p(np.abs(database.pressure), database.latitude) # compute depth from pressure

Zonal averages of surface tritium

Read model surface data

In [5]:
ds_model = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/surface.%s.nc'%case, decode_times=False)
ds_model = fixtime(ds_model, 1824)
ds_model = ds_model.squeeze(dim='z_t')
ds_model = ds_model.drop(['z_t'])
ds_model['TAREA'] = ds_model.TAREA.where(ds_model.REGION_MASK>0)     # mask land & marginal seas
ds_model['TRITIUM'] = ds_model.TRITIUM.where(ds_model.REGION_MASK>0)
ds_model['TRITIUM'] = xr.where(ds_model.TRITIUM<0, 0, ds_model.TRITIUM) # remove negative values
ds_model['time'] = ds_model.time.to_index().to_datetimeindex()

Surface tritium by latitude bands

Shaded area corresponds to 95% confidence interval for model mean.

In [6]:
fig, axs = plt.subplots(len(lat_bands),1,sharex=True,figsize=(6.4,10))
fig.subplots_adjust(top=0.95,hspace=0.3)
for ax, band in zip(axs.ravel(), lat_bands[::-1]):
    df_obs = get_lat_band(database, *band)
    df_obs['tritium'] = df_obs.tritium * 0.10995 # TU -> pmol/m^3
    df_obs.rename({'tritium': 'observations'}, axis='columns', inplace=True)
    df_mod = lat_band_mean(ds_model,*band)
    _ = df_obs.plot(x='stadate', y='observations', style='.', ax=ax, legend=False, markersize=2)
    _ = df_mod['model'].plot(ax=ax, style='-', legend=False)
    _ = ax.fill_between(df_mod.index, df_mod.vmin, df_mod.vmax, facecolor='C1', alpha=0.3)
    _ = ax.set_title('lat: {} to {}'.format(*band), fontsize=10)
    _ = ax.legend(loc='best')
    ax.autoscale(axis='x',tight=True)
    
_ = fig.text(0.05,0.575,'surface tritium (pmol m$^{-3}$)',ha='center',va='center',rotation=90)

Zonal averages of tritium inventory

Read model tritium inventory

In [7]:
ds_model_inv = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/tritium.%s.nc'%case)
ds_model_inv['TAREA'] = ds_model_inv.TAREA.where(ds_model_inv.REGION_MASK>0) # mask land & marginal seas
ds_model_inv['time'] = ds_model_inv.time.to_index().to_datetimeindex()

Tritium inventory by latitude bands

Shaded area corresponds to 95% confidence interval for model mean.

In [8]:
fig, axs = plt.subplots(len(lat_bands),1,sharex=True,figsize=(6.4,10))
fig.subplots_adjust(top=0.95,hspace=0.3)
for ax, band in zip(axs.ravel(), lat_bands[::-1]):
    df_obs = get_lat_band(database, *band, zsurf=None)
    profile, x, y = get_profile(df_obs, zb=zbins2)
    inventory = get_inventory(profile,nn=11,zb=zbins2)
    df_mod = lat_band_mean(ds_model_inv,*band,'TRITIUM_inventory')
    _ = inventory.plot(style='.', ax=ax, legend=False, markersize=2)
    _ = df_mod['model'].plot(ax=ax, style='-', legend=False)
    _ = ax.fill_between(df_mod.index, df_mod.vmin, df_mod.vmax, facecolor='C1', alpha=0.3)
    _ = ax.set_title('lat: {} to {}'.format(*band), fontsize=10)
    _ = ax.legend(loc='best')
    ax.autoscale(axis='x',tight=True)
    
_ = fig.text(0.05,0.575,'tritium inventory (pmol m$^{-2}$)',ha='center',va='center',rotation=90)

Model zonal averages of tritium inventory

In [9]:
from scipy.interpolate import griddata
from ccsm_utils import grid_area
from ipyparallel import Client

cl = Client()
nprocs = len(cl.ids)
if nprocs > 0: print('running on %d cores'%nprocs)
dview = cl[:]
dview.block = True
dview.execute('import numpy as np')
dview.execute('from scipy.interpolate import griddata')

lview = cl.load_balanced_view()
lview.block = True

def pregrid(data):
    work1 = np.array(data.tolist(),copy=False)
    work2 = griddata((lonold,latold),work1,(lonnew,latnew),method='linear')
    work3 = np.ma.masked_where(np.isnan(work2),work2)
    return np.ma.average(work3, axis=1, weights=area)

running on 24 cores

Atlantic

In [10]:
lon, lat = np.arange(-97,19,1), np.arange(-33.5,69,0.5)
ulon, ulat = np.hstack([lon - 0.5, lon[-1] + 0.5]), np.hstack([lat - 0.25, lat[-1] + 0.25])
lon_new, lat_new = np.meshgrid(lon, lat)
area = grid_area(ulon,ulat)

tlon = xr.where(ds_model_inv.TLONG>180,ds_model_inv.TLONG-360,ds_model_inv.TLONG).values
tlat = ds_model_inv.TLAT.values
data = ds_model_inv.TRITIUM_inventory.where(ds_model_inv.REGION_MASK==6).values

dview.push(dict(lonold=tlon.ravel(),latold=tlat.ravel(),lonnew=lon_new,latnew=lat_new,area=area))

ntime = ds_model_inv.TRITIUM_inventory.shape[0]
arglist = [data[t,...].ravel() for t in range(ntime)]
result = lview.map_sync(pregrid,arglist)
zavg_atl = np.ma.array(result)
In [11]:
# nlon, nlat = lon.shape[0], lat.shape[0]
# zavg_atl = np.ma.masked_all((ntime,nlat),dtype=np.float)

# for n in range(ntime):
#     work = griddata((tlon.data.ravel(),tlat.data.ravel()),data[n,...].data.ravel(),(lon_new,lat_new),method='linear')
#     work = np.ma.masked_where(np.isnan(work),work)
#     zavg_atl[n,...] = np.ma.average(work, axis=1, weights=area)
In [12]:
fig, ax = plt.subplots()
pm = ax.pcolormesh(ds_model_inv.time, lat, zavg_atl.transpose(),cmap=plt.cm.plasma)
_ = ax.set(xlabel='time', ylabel='latitude', title='Atlantic tritium inventory zonal average')
cb = fig.colorbar(pm,ax=ax)
cb.set_label('pmol m$^{-2}$')

Pacific

In [13]:
lon, lat = np.arange(100,291,1), np.arange(-33.5,67,0.5)
ulon, ulat = np.hstack([lon - 0.5, lon[-1] + 0.5]), np.hstack([lat - 0.25, lat[-1] + 0.25])
lon_new, lat_new = np.meshgrid(lon, lat)
area = grid_area(ulon,ulat)

tlon = ds_model_inv.TLONG.values
tlat = ds_model_inv.TLAT.values
data = ds_model_inv.TRITIUM_inventory.where(ds_model_inv.REGION_MASK==2).values

dview.push(dict(lonold=tlon.ravel(),latold=tlat.ravel(),lonnew=lon_new,latnew=lat_new,area=area))

ntime = ds_model_inv.TRITIUM_inventory.shape[0]
arglist = [data[t,...].ravel() for t in range(ntime)]
result = lview.map_sync(pregrid,arglist)
zavg_pac = np.ma.array(result)
In [14]:
# nlon, nlat = lon.shape[0], lat.shape[0]
# zavg_pac = np.ma.masked_all((ntime,nlat),dtype=np.float)

# for n in range(ntime):
#     work = griddata((tlon.data.ravel(),tlat.data.ravel()),data[n,...].data.ravel(),(lon_new,lat_new),method='linear')
#     work = np.ma.masked_where(np.isnan(work),work)
#     zavg_pac[n,...] = np.ma.average(work, axis=1, weights=area)
In [15]:
fig, ax = plt.subplots()
pm = ax.pcolormesh(ds_model_inv.time, lat, zavg_pac.transpose(),cmap=plt.cm.plasma)
_ = ax.set(xlabel='time', ylabel='latitude', title='Pacific tritium inventory zonal average')
cb = fig.colorbar(pm,ax=ax)
cb.set_label('pmol m$^{-2}$')