In [1]:
import warnings
warnings.filterwarnings('ignore')
import numpy
import pandas
import iris
import iris.coord_categorisation
import matplotlib.pyplot as plt
import glob
import scipy
In [2]:
%matplotlib inline
In [3]:
#tfile = '/g/data/ua6/DRSv3/CMIP5/CCSM4/historicalGHG/mon/ocean/r1i1p1/thetao/latest/thetao_Omon_CCSM4_historicalGHG_r1i1p1_200001-200512.nc'
#vfile = '/g/data/r87/dbi599/DRSv2/CMIP5/CCSM4/historical/fx/ocean/r0i0p0/volcello/latest/volcello-inferred_fx_CCSM4_historical_r0i0p0.nc'
#bfile = '/g/data/r87/dbi599/DRSv2/CMIP5/CCSM4/historical/fx/ocean/r0i0p0/basin/latest/basin_fx_CCSM4_historical_r0i0p0.nc'
tfile = '/g/data/r87/dbi599/data_argo/raw/to-anom_Oyr_ArgoRoemmich_2004-2015.nc'
tclim_file = '/g/data/r87/dbi599/data_argo/raw/to-annual-clim_Omon_ArgoRoemmich_2004-2015.nc'
vfile = '/g/data/r87/dbi599/data_argo/raw/volcello-inferred_fx_ArgoRoemmich.nc'
In [4]:
level_subset = lambda cell: cell <= 2000
level_constraint = iris.Constraint() #depth=level_subset)
tanom_cube = iris.load_cube(tfile, level_constraint)
tclim_cube = iris.load_cube(tclim_file)
vcube = iris.load_cube(vfile, level_constraint)
#bcube = iris.load_cube(bfile ,level_constraint)
In [5]:
tanom_cube
Out[5]:
In [6]:
tclim_cube
Out[6]:
In [7]:
#tcube.collapsed('time', iris.analysis.MEAN)
#tcube.data = tcube.data - 273.15
tcube = tclim_cube.copy()
tcube.data = tanom_cube.collapsed('time', iris.analysis.MEAN).data + tclim_cube.data
tcube
Out[7]:
In [8]:
tdata = tcube.data.flatten()
vdata = vcube.data.flatten()
#bdata = bcube.data.flatten()
In [9]:
df = pandas.DataFrame(index=range(tdata.shape[0]))
df['temperature'] = tdata.filled(fill_value=5000)
df['volume'] = vdata.filled(fill_value=5000)
#df['basin'] = bdata.filled(fill_value=5000)
In [10]:
df = df[df.temperature != 5000]
#df = df[df.temperature != -273.15]
In [11]:
print(df['temperature'].values.min())
print(df['temperature'].values.max())
In [12]:
bin_edges = numpy.arange(-2.5, 31.5, 1)
bin_edges
Out[12]:
In [13]:
dvdt, edges, binnum = scipy.stats.binned_statistic(df['temperature'].values, df['volume'].values, statistic='sum', bins=bin_edges)
In [14]:
temps = (edges[1:] + edges[:-1]) / 2
temps
Out[14]:
In [15]:
plt.figure(figsize=(12,7))
plt.plot(temps, dvdt)
plt.xlabel('temperature (C)')
plt.ylabel('$m^3 / C$')
plt.title('Volumetic distribution (dV/dT)')
plt.show()
In [16]:
vt = dvdt.cumsum()
In [17]:
plt.plot(temps, vt)
plt.xlabel('temperature (C)')
plt.ylabel('$m^3 / C$')
plt.title('V(T)')
plt.show()
V(T) is the volume of water colder than T.
In [18]:
#infiles = glob.glob('/g/data/ua6/DRSv3/CMIP5/CCSM4/historicalGHG/mon/ocean/r1i1p1/thetao/latest/thetao_Omon_CCSM4_historicalGHG_r1i1p1_??????-??????.nc')
#infiles.sort()
In [19]:
tanom_cube
Out[19]:
In [21]:
vt_timeseries = numpy.array([])
#for infile in infiles:
# print(infile)
# tcube = iris.load_cube(infile, level_constraint)
#iris.coord_categorisation.add_year(tanom_cube, 'time')
#iris.coord_categorisation.add_month(tanom_cube, 'time')
#tcube = tanom_cube.aggregated_by(['year'], iris.analysis.MEAN)
#tcube.remove_coord('year')
#tcube.remove_coord('month')
for time_slice in tanom_cube.slices_over('time'):
tdata = time_slice.data + tclim_cube.data
tdata = tdata.flatten()
#tdata = tdata - 273.15
df = pandas.DataFrame(index=range(tdata.shape[0]))
df['temperature'] = tdata.filled(fill_value=5000)
df['volume'] = vdata.filled(fill_value=5000)
df = df[df.temperature != 5000]
#df = df[df.temperature != -273.15]
dvdt, edges, binnum = scipy.stats.binned_statistic(df['temperature'].values, df['volume'].values, statistic='sum', bins=bin_edges)
vt = dvdt.cumsum()
vt_timeseries = numpy.vstack([vt_timeseries, vt]) if vt_timeseries.size else vt
In [22]:
vt_timeseries.shape
Out[22]:
In [23]:
def linear_trend(data, time_axis):
"""Calculate the linear trend.
polyfit returns [b, a] corresponding to y = a + bx
"""
return numpy.polyfit(time_axis, data, 1)[0]
ntime = vt_timeseries.shape[0]
seconds = numpy.arange(ntime) * 60 * 60 * 24 * 365.25
trends = numpy.apply_along_axis(linear_trend, 0, vt_timeseries, seconds)
In [24]:
plt.plot(temps, trends)
plt.xlabel('temperature (C)')
plt.ylabel('$m^3 / s$')
plt.title('Transformation rate (dV/dt)')
plt.show()
It makes sense to me that the transformation rates would be negative. In a warming climate, the volume of water cooler than a given temperature would likely decrease. This doesn't match Jan's Argo results, which shows positive rates...
In [26]:
plt.figure(figsize=(12,7))
plt.plot(temps, trends / dvdt)
plt.xlabel('temperature (C)')
plt.ylabel('$C / s$')
plt.title('Implied diabatic temperature tendency (dV/dt / dV/dT)')
plt.show()
In [ ]: