In [1]:
__depends__ = ['../outputs/llc_kuroshio_timeseries.nc','../outputs/VelStats_anomaly_24h_4320_0_spectral.npz',
'../outputs/VelStats_anomaly_24h_spectral.npz']
__dest__ = []
This notebook showcases the analysis applied to LLC outputs. Here the calculations are performed for a single snapshot. The full LLC model outputs can be obtained from the ECCO Project. All fields used in this paper take about 700 GB!
The analysis leverage on other pieces of code developed by the first author: llctools and pyspec.
In [2]:
import datetime
import numpy as np
import scipy as sp
from scipy import interpolate
import matplotlib.pyplot as plt
%matplotlib inline
import cmocean
import seawater as sw
from netCDF4 import Dataset
from llctools import llc_model
from pyspec import spectrum as spec
In [3]:
c1 = 'slateblue'
c2 = 'tomato'
c3 = 'k'
c4 = 'indigo'
plt.rcParams['lines.linewidth'] = 1.5
ap = .75
plt.style.use('seaborn-colorblind')
def leg_width(lg,fs):
"""" Sets the linewidth of each legend object """
for legobj in lg.legendHandles:
legobj.set_linewidth(fs)
def parse_time(times):
"""
Converts an array of strings that defines
the LLC outputs into datatime arrays,
e.g., '20110306T010000' --> datetime.datetime(2011, 3, 6, 1, 0)
Input
------
times: array of strings that define LLC model time
Output
------
time: array of datetime associated with times
"""
time = []
for i in range(times.size):
yr = times[i][:4]
mo = times[i][4:6]
day = times[i][6:8]
hr = times[i][9:11]
time.append(datetime.datetime(int(yr),int(mo),int(day),int(hr)))
return np.array(time)
In [4]:
grid_path = '../data/llc/2160/grid/'
data_path = '../data/llc/2160/uv/'
In [5]:
# Kuroshio Extension model class
m = llc_model.LLCRegion(grid_dir=grid_path,data_dir=data_path,Nlon=480,Nlat=466,Nz=1)
In [6]:
m.load_grid()
In [7]:
# model sub-region surface fields files
fileu = m.data_dir+'U_480x466x1.20110308T220000'
filev = m.data_dir+'V_480x466x1.20110308T220000'
fileeta = m.data_dir[:-3]+'Eta/Eta_480x466x1.20110308T220000'
In [8]:
time_string = fileu[-15:]
time=llc_model.parse_time(time_string)
time
Out[8]:
In [9]:
# important note: U,V are relative to the LLC model grid,
# not geostrophical coordinates. Thus, on
# faces 4 and 5, U = meridional component
# and V = -zonal component (see Dimitris's llc.readme).
u, v, eta = m.load_2d_data(filev), -m.load_2d_data(fileu), m.load_2d_data(fileeta)
In [10]:
lon,lat = m.lon[m.Nlat//2],m.lat[:,m.Nlon//2]
# create a regular Cartesian grid
dd = 6. # grid spacing [km]
dlon = dd/111.320*np.cos(np.abs(m.lat[m.Nlat//2,m.Nlon//2])*np.pi/180.)
dlat = dd/110.574
lonimin,lonimax = lon.min()+dlon,lon.max()-dlon
latimin,latimax = lat.min()+dlat,lat.max()-dlat
loni = np.arange(m.lon.min(),m.lon.max()+dlon,dlon)
lati = np.arange(m.lat.min(),m.lat.max()+dlat,dlat)
long,latg = np.meshgrid(loni,lati)
f0 = sw.f(latg)
In [11]:
interpu, interpv, interpeta = sp.interpolate.interp2d(lon,lat,u), sp.interpolate.interp2d(lon,lat,v), sp.interpolate.interp2d(lon,lat,eta)
In [12]:
ui, vi,etai = interpu(loni,lati), interpv(loni,lati), interpeta(loni,lati)
In [13]:
def calc_gradu(u,v,dd = 6.):
uy,ux = np.gradient(u,dd,dd)
vy,vx = np.gradient(v,dd,dd)
vort, div, strain = (vx - uy), ux+vy, ( (ux-vy)**2 + (vx+uy)**2 )**.5
return vort, div, strain
# double mirror ui and vi
def double_mirror(a,forward='True'):
if forward:
A = np.hstack([a,np.fliplr(a)])
A = np.vstack([A,np.fliplr(A)])
else:
iy,ix = a.shape
A = a[:iy//2,:ix//2]
return A
def calc_gradu2(u,v,dd = 6.):
u, v = double_mirror(u), double_mirror(v)
iy,ix = u.shape
Lx, Ly = (ix-1)*dd, (iy-1)*dd
dk = 1./Lx
dl = 1./Ly
l = 2*np.pi*dl*np.append( np.arange(0.,iy//2), np.arange(-iy//2,0.) )
k = 2*np.pi*dk*np.arange(0.,ix//2+1)
k,l = np.meshgrid(k,l)
uh, vh = np.fft.rfft2(u), np.fft.rfft2(v)
ux, uy = np.fft.irfft2(1j*k*uh), np.fft.irfft2(1j*l*uh)
vx, vy = np.fft.irfft2(1j*k*vh), np.fft.irfft2(1j*l*vh)
vort, div, strain = (vx - uy), ux+vy, ( (ux-vy)**2 + (vx+uy)**2 )**.5
return vort, div, strain
def rms(field):
return ((field**2).mean())**.5
In [14]:
vort, div, strain = calc_gradu(ui,vi,dd = 6.e3)
vort, div, strain = vort/f0, div/f0, strain/f0
In [15]:
vort2, div2, strain2 = calc_gradu2(ui,vi,dd = 6.e3)
vort2,div2, strain2 = double_mirror(vort2,forward=False),double_mirror(div2,forward=False), double_mirror(strain2,forward=False)
vort2, div2, strain2 = vort2/f0, div2/f0, strain2/f0
In [16]:
vort.mean()/np.abs(vort).max(), div.mean()/np.abs(div).max(), strain.mean()/np.abs(strain).max()
Out[16]:
In [17]:
vort2.mean()/np.abs(vort2).max(), div2.mean()/np.abs(div2).max(), strain2.mean()/np.abs(strain2).max()
Out[17]:
In [18]:
fig = plt.figure(figsize=(14,4))
cv = np.linspace(-1.5,1.5,20)
cd = np.linspace(-.5,.5,20)
cs = np.linspace(0.,1.5,10)
ax = fig.add_subplot(131)
plt.contourf(vort,cv,vmin=cv.min(),vmax=cv.max(),cmap=cmocean.cm.balance,extend='both')
plt.title('vorticity, rms = %f' % rms(vort))
#plt.colorbar()
plt.xticks([]); plt.yticks([])
ax = fig.add_subplot(132)
plt.contourf(vort2,cv,vmin=cv.min(),vmax=cv.max(),cmap=cmocean.cm.balance,extend='both')
plt.title('vorticity, rms = %f' % rms(vort2))
#plt.colorbar()
plt.xticks([]); plt.yticks([])
Out[18]:
In [19]:
fig = plt.figure(figsize=(14,4))
ax = fig.add_subplot(131)
plt.contourf(div,cd,vmin=cd.min(),vmax=cd.max(),cmap=cmocean.cm.balance,extend='both')
plt.title('divergence, rms = %f' % rms(div))
#plt.colorbar()
plt.xticks([]); plt.yticks([])
ax = fig.add_subplot(132)
plt.contourf(div2,cd,vmin=cd.min(),vmax=cd.max(),cmap=cmocean.cm.balance,extend='both')
plt.title('divergence, rms = %f' % rms(div2))
#plt.colorbar()
plt.xticks([]); plt.yticks([])
Out[19]:
In [20]:
fig = plt.figure(figsize=(14,4))
ax = fig.add_subplot(131)
plt.contourf(strain,cs,vmin=cs.min(),vmax=cs.max(),cmap=cmocean.cm.amp,extend='both')
plt.title('divergence, rms = %f' % rms(strain))
#plt.colorbar()
plt.xticks([]); plt.yticks([])
ax = fig.add_subplot(132)
plt.contourf(strain2,cs,vmin=cs.min(),vmax=cs.max(),cmap=cmocean.cm.amp,extend='both')
plt.title('strain, rms = %f' % rms(strain2))
#plt.colorbar()
plt.xticks([]); plt.yticks([])
Out[20]:
In [21]:
stats_4320 = np.load(__depends__[1])
stats_2160 = np.load(__depends__[2])
llc = Dataset(__depends__[0])
In [22]:
time2160 = parse_time(llc['2160']['hourly']['time'][:])
timed2160 = time2160[::24]
time4320 = parse_time(llc['4320']['hourly']['time'][:])
timed4320 = time4320[::24]
In [23]:
cv = np.linspace(-1.5,1.5,20)
cd = np.linspace(-.5,.5,20)
cs = np.linspace(0.,1.5,10)
fig = plt.figure(figsize=(19,4))
ax = fig.add_subplot(131)
plt.contourf(vort,cv,vmin=cv.min(),vmax=cv.max(),cmap='RdBu_r',extend='both')
plt.title('vorticity, rms = %f' % rms(vort))
plt.colorbar()
plt.xticks([]); plt.yticks([])
ax = fig.add_subplot(132)
plt.title('divergence, rms = %f' % rms(div))
plt.contourf(div,cd,vmin=cd.min(),vmax=cd.max(),cmap='RdBu_r',extend='both')
plt.colorbar()
plt.xticks([]); plt.yticks([])
ax = fig.add_subplot(133)
plt.title('strain rate, rms %f' % rms(strain))
plt.contourf(strain,cs,vmax=cs.max(),cmap='viridis',extend='max')
plt.colorbar()
plt.xticks([]); plt.yticks([])
Out[23]:
In [24]:
specU = spec.TWODimensional_spec(ui.copy(),d1=dd,d2=dd)
specV = spec.TWODimensional_spec(vi.copy(),d1=dd,d2=dd)
specEta = spec.TWODimensional_spec(etai.copy(),d1=dd,d2=dd)
In [25]:
iEu,iEv, iEeta = specU.ispec,specV.ispec, specEta.ispec
iE = 0.5*(iEu+iEv)
In [26]:
kr = np.array([1.e-4,1.])
e2 = kr**-2/1.e4
e3 = kr**-3/1.e7
e5 = kr**-5/1.e9
fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(121)
plt.loglog(specU.ki,iE)
plt.loglog(kr,12.*e2,'.5',linewidth=2); plt.text(1/17.5,5.e-1,'-2',fontsize=14)
plt.loglog(kr,35*e3,'.5',linewidth=2); plt.text(1/30.,2.e-2,'-3',fontsize=14)
plt.xlim(1.e-3,1.e-1)
plt.ylim(1.e-2,1.e2)
plt.xlabel('Wavenumber [cpkm]')
plt.ylabel(r'KE density [m$^2$ s$^{-2}$/cpkm]')
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=.45, hspace=None)
ax = fig.add_subplot(122)
plt.loglog(specEta.ki,iEeta)
plt.loglog(kr,e2/.5e1,'.5',linewidth=2); plt.text(1/17.5,1.e-2,'-2',fontsize=14)
plt.loglog(kr,3*e5/1.5e2,'.5',linewidth=2); plt.text(1/25.5,1.e-5,'-5',fontsize=14)
plt.xlim(1.e-3,1.e-1)
plt.ylim(1.e-6,1.e2)
plt.ylabel(r'SSH variance density [m$^2$/cpkm]')
plt.xlabel('Wavenumber [cpkm]')
Out[26]:
In [27]: