Compute column inventories of tritium at HOTS

Ivan Lima - Wed Apr 11 2018 14:52:41 -0400

%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 bgc_utils import rho
from scipy import interpolate
from he3_inventory import *
pd.options.display.max_columns = 50

Read database

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

Extract observational data at HOTS location & compute tritium inventory

hots_lat = database[database.sect_id=='hot'].latitude.mean()
hots_lon = database[database.sect_id=='hot'].longitude.mean()
hots = get_box(database, hots_lon, hots_lat, dlon=6, dlat=3) # all data in a 6x3 degree box

zbins2 = np.array([-4000, -3500, -3000, -2500, -2000, -1500, -1000, -500, -400, -300, -250, -200, -150, -100, -50, 0])
hots_profile, x, y = get_profile(hots, zb=zbins2)
hots_inventory = get_inventory(hots_profile,nn=11,zb=zbins2)

Read model output at HOTS & compute tritium inventory

ds_hots = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/HOT.tr3he.GIAF.01.nc', decode_times=False)
ds_hots = fixtime(ds_hots, 1824)
ds_hots = ds_hots.squeeze(dim=('nlat','nlon'))
ds_hots = ds_hots.drop(['TLONG','TLAT'])
model_tritium_inv = (ds_hots.TRITIUM * ds_hots.dz * 1.e-2).sum(dim='z_t').to_series() # cm -> m

Plot figures

Tritium profiles & inventory from observations

fig, axs = plt.subplots(2,1,sharex=True,figsize=(9.5,9))
fig.subplots_adjust(hspace=0.05,right=0.875)
_ = axs[0].plot(hots_inventory.index, hots_inventory.values, '.', label='observations')
_ = axs[0].set(ylabel='tritium inventory (pmol m$^{-2}$)', title='HOTS')

pm = axs[1].pcolormesh(x, y, hots_profile, cmap=plt.cm.plasma, norm=colors.LogNorm())
_ = axs[1].set(xlabel='time', ylabel='depth (m)')
l, b, w, h =  axs[1].get_position().bounds
cax = fig.add_axes((l+w+0.02, b, 0.02, h))
cb = fig.colorbar(pm, cax=cax)
cb.set_label('tritium concentration (pmol m$^{-3}$)')
#_ = axs[1].plot(hots.stadate, hots.depth, 'x')

Tritium inventory from model and observations

fig, ax = plt.subplots()
_ = hots_inventory.plot(ax=ax, style='.', label='observations')
# _ = hots_inventory.resample('A').mean().plot(ax=ax, style='.', label='observations')
_ = model_tritium_inv.plot(ax=ax, style='.', label='model')
_ = ax.set(ylabel='tritium inventory (pmol m$^{-2}$)', title='HOTS')
_ = ax.legend(loc='best')

Surface tritium from model and observations

hots_surface = hots_profile.iloc[-1,:]
model_tritium_surf = ds_hots.TRITIUM.sel(z_t=500).to_series()
fig, ax = plt.subplots()
_ = hots_surface.plot(ax=ax, style='.', label='observations')
_ = model_tritium_surf.plot(ax=ax, style='.', label='model')
_ = ax.set(ylabel='surface tritium (pmol m$^{-3}$)', title='HOTS')
_ = ax.legend(loc='best')

Tritium profiles from model and observations

fig, axs = plt.subplots(2,1,sharex=True,figsize=(9.5,9))
fig.subplots_adjust(right=0.875)
pm = axs[0].pcolormesh(x, y, hots_profile, cmap=plt.cm.plasma, norm=colors.LogNorm())
_ = axs[0].set(ylabel='depth (m)', title='observations')
l, b, w, h =  axs[0].get_position().bounds
cax = fig.add_axes((l+w+0.02, b, 0.02, h))
cb = fig.colorbar(pm, cax=cax)
cb.set_label('tritium (pmol m$^{-3}$)')
ymin, ymax = axs[0].get_ylim()
vmin, vmax = pm.get_clim()
data = da2ma(ds_hots.TRITIUM.transpose().values)
data[data<=0] = vmin
pm = axs[1].pcolormesh(ds_hots.time, -ds_hots.z_t/100, data, cmap=plt.cm.plasma, vmin=vmin, vmax=vmax, 
                       norm=colors.LogNorm())
_ = axs[1].set(xlabel='time', ylabel='depth (m)', ylim=(ymin,ymax), title='model')