import numpy as np
from scipy import io
import matplotlib.pyplot as plt
%matplotlib inline
Vendor: Continuum Analytics, Inc. Package: mkl Message: trial mode expires in 13 days Vendor: Continuum Analytics, Inc. Package: mkl Message: trial mode expires in 13 days
plt.rcParams.update({'font.size': 12
, 'legend.markerscale': 1., 'axes.titlesize': 12, 'axes.labelsize' : 12,
'legend.fontsize' : 10,'legend.handlelength': 3})
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
color2 = '#6495ed'
color1 = '#ff6347'
color5 = '#8470ff'
color3 = '#3cb371'
color4 = '#ffd700'
color6 = '#ba55d3'
lw1=3
aph=.7
data_path = './outputs/'
spec = np.load(data_path+"isotropic_ke_spectra_0m.npz")
spec1000 = np.load(data_path+"isotropic_ke_spectra_1000m.npz")
ssh = np.load(data_path+"isotropic_spectrum_eta.npz")
## -2 and -3 slopes in the loglog space
ks = np.array([1.e-3,1])
Es2 = .2e-4*(ks**(-2))
Es3 = .5e-6*(ks**(-3))
rd1 = 22.64 # [km]
Enoise = np.ones(2)*2.*1.e-4
def add_second_axis(ax1):
""" Add a x-axis at the top of the spectra figures """
ax2 = ax1.twiny()
ax2.set_xscale('log')
ax2.set_xlim(ax1.axis()[0], ax1.axis()[1])
kp = 1./np.array([500.,200.,100.,40.,20.,10.,5.])
lp=np.array([500,200,100,40,20,10,5])
ax2.set_xticks(kp)
ax2.set_xticklabels(lp)
plt.xlabel('Wavelength [km]')
Surface spectrum
lw,aph=3,.36
fig = plt.figure(facecolor='w', figsize=(12.,10.))
ax1 = fig.add_subplot(111)
ax1.fill_between(spec['k'],spec['El'],spec['Eu'],\
color=color1, alpha=aph)
ax1.fill_between(spec['k'],spec['Efl'],spec['Efu'],\
color=color1, alpha=aph)
ax1.fill_between(spec['k'],spec['Endl'],spec['Endu'],\
color=color3, alpha=aph)
ax1.fill_between(spec['k'],spec['Edl'],spec['Edu'],\
color=color4, alpha=aph)
ax1.fill_between(spec['k'],spec['Endfl'],spec['Endfu'],\
color=color3, alpha=aph)
ax1.fill_between(spec['k'],spec['Edfl'],spec['Edfu'],\
color=color4, alpha=aph)
ax1.fill_between(spec['kg'],spec['Egl'],spec['Egu'],\
color=color2, alpha=aph)
ax1.fill_between(spec['kg'],spec['Egfl'],spec['Egfu'],\
color=color2, alpha=aph)
ax1.loglog(spec['k'],spec['E'],color=color1,linewidth=lw,
label=r'total')
ax1.loglog(spec['k'],spec['Ef'],'--',color=color1,linewidth=lw)
ax1.loglog(spec['k'],spec['End'],color=color3,linewidth=lw,
label=r'rotational')
ax1.loglog(spec['k'],spec['Endf'],'--',color=color3,linewidth=lw)
ax1.loglog(spec['k'],spec['Ed'],color=color4,linewidth=lw,
label=r'divergent')
ax1.loglog(spec['k'],spec['Edf'],'--',color=color4,linewidth=lw)
ax1.loglog(spec['kg'],spec['Eg'],color=color2,linewidth=lw,
label=r'from SSH')
ax1.loglog(spec['kg'],spec['Egf'],'--',color=color2,linewidth=lw)
ax1.loglog(ks,Es2,'--', color='k',linewidth=2.,alpha=.7)
ax1.loglog(ks,Es3,'--', color='k',linewidth=2.,alpha=.7)
plt.text(0.0025, 1.941,u'k$^{-2}$')
plt.text(0.0060, 1.951,u'k$^{-3}$')
plt.xlabel('Wavenumber [cpkm]')
plt.ylabel(u'KE spectral density [ m$^{2}$ s$^{-2}$/ cpkm]')
plt.text(1./5, 3., "d", size=35, rotation=0.)
plt.text(1./20., 3., "llc 4320 simulation", size=25, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
lg = plt.legend(loc=3,title='isotropic spectrum, 0 m', numpoints=1,ncol=2)
lg.draw_frame(False)
plt.axis((1./1.e3,1./4.,.5/1.e4,5e0))
add_second_axis(ax1)
plt.savefig('figs/iso_spec_model_0m',bbox_inches='tight')
lw,aph=3,.36
fig = plt.figure(facecolor='w', figsize=(12.,10.))
ax1 = fig.add_subplot(111)
ax1.fill_between(spec1000['k'],spec1000['El'],spec1000['Eu'],\
color=color1, alpha=aph)
ax1.fill_between(spec1000['k'],spec1000['Efl'],spec1000['Efu'],\
color=color1, alpha=aph)
ax1.fill_between(spec1000['k'],spec1000['Endl'],spec1000['Endu'],\
color=color3, alpha=aph)
ax1.fill_between(spec1000['k'],spec1000['Edl'],spec1000['Edu'],\
color=color4, alpha=aph)
ax1.fill_between(spec1000['k'],spec1000['Endfl'],spec1000['Endfu'],\
color=color3, alpha=aph)
ax1.fill_between(spec1000['k'],spec1000['Edfl'],spec1000['Edfu'],\
color=color4, alpha=aph)
ax1.loglog(spec1000['k'],spec1000['E'],color=color1,linewidth=lw,
label=r'total')
ax1.loglog(spec1000['k'],spec1000['Ef'],'--',color=color1,linewidth=lw)
ax1.loglog(spec1000['k'],spec1000['End'],color=color3,linewidth=lw,
label=r'rotational')
ax1.loglog(spec1000['k'],spec1000['Endf'],'--',color=color3,linewidth=lw)
ax1.loglog(spec1000['k'],spec1000['Ed'],color=color4,linewidth=lw,
label=r'divergent')
ax1.loglog(spec1000['k'],spec1000['Edf'],'--',color=color4,linewidth=lw)
ax1.loglog(ks,Es2,'--', color='k',linewidth=2.,alpha=.7)
ax1.loglog(ks,Es3,'--', color='k',linewidth=2.,alpha=.7)
plt.text(0.0025, 1.941,u'k$^{-2}$')
plt.text(0.0060, 1.951,u'k$^{-3}$')
plt.xlabel('Wavenumber [cpkm]')
plt.ylabel(u'KE spectral density [ m$^{2}$ s$^{-2}$/ cpkm]')
plt.text(1./5, 3., "e", size=35, rotation=0.)
plt.text(1./20., 3., "llc 4320 simulation", size=25, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
lg = plt.legend(loc=3,title='isotropic spectrum, 1000 m', numpoints=1,ncol=2)
lg.draw_frame(False)
plt.axis((1./1.e3,1./4.,.5/1.e4,5e0))
add_second_axis(ax1)
plt.savefig('figs/iso_spec_model_1000m',bbox_inches='tight')
The SSH variance spectrum of hourly fields presents a "break" at about 40 km; the shape of the spectrum changes from about $\kappa^{-5}$ to $\kappa^{-1}$. We attribute this signigicant spectral flattening to the contribution of unbalanced flows, likely dominated by inertia-gravity waves. These motions are very incoherent. Removing tidal frequencies using a python implementation of T_TIDE does not affect the spectrum at small scales; tidal motions are dominated by barotropic tides and project onto very large scales.
Filtering out all fast flows by daily-averaging supress most of the high-wavenumber variance. The spectrum of filtered motions is consistent with predictions of isotropic interior QG turbulence.
import seaborn as sns
#sns.set(style="darkgrid")
sns.set(style="whitegrid")
from pyspec import spectrum
from plots_aux import *
/Users/crocha/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter. warnings.warn(self.msg_depr % (key, alt_key)) Vendor: Continuum Analytics, Inc. Package: mkl Message: trial mode expires in 13 days
lw = 1.
aph=.35
sn = 11
L = 1./ssh['k']
f1 = L > 100.
f2 = (L<=100.)&(L>20.)
f3 = L<=20.
SN = np.zeros_like(ssh['k'])
SN[f1],SN[f2],SN[f3] = sn, 4*sn,8*sn
Y = np.ones_like(SN)*.25e-4
win = np.hanning(spec['k'].size)
fac = (spec['k'].size/(win**2).sum())
fig = plt.figure(figsize=(8.27/2-.25,11.69/3-.25))
ax1 = fig.add_subplot(111)
ax1.fill_between(ssh['k'],fac*ssh['El'],fac*ssh['Eu'],\
color=color1, alpha=aph)
ax1.fill_between(ssh['k'],fac*ssh['Efl'],fac*ssh['Efu'],\
color=color2, alpha=aph)
ax1.fill_between(ssh['k'],fac*ssh['Edt_l'],fac*ssh['Edt_u'],\
color=color3, alpha=aph)
ax1.loglog(ssh['k'],fac*ssh['E'],color=color1,linewidth=lw,
label=r'hourly')
ax1.loglog(ssh['k'],fac*ssh['Ef'],color=color2,linewidth=lw,
label=r'daily-averaged')
ax1.loglog(ssh['k'],fac*ssh['Edt'],color=color3,linewidth=lw,
label=r'de-tided')
ax1.loglog(ks,Es2,'-', color='0.5',linewidth=1.)
ax1.loglog(ks,Es3,'-', color='0.5',linewidth=1.)
plt.text(0.0013, 4.25,u'-2')
plt.text(0.00475, 4.25,u'-3')
plt.xlabel('Wavenumber [cpkm]')
plt.ylabel(u'SSH variance density [m$^{2}$ / cpkm]')
#plt_spec_error(x=ssh['k'],y=Y,sn=SN)
lg = plt.legend(loc=3,title='', numpoints=1,ncol=1)
lg.draw_frame(False)
plt.text(1./20., 5, "llc 4320 simulation", size=10, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
plt.axis((1./1.e3,1./4.,1./1.e7,1.6e1))
add_second_axis(ax1)
plt.savefig('figs/iso_spec_model_ssh',bbox_inches='tight')
plt.savefig('figs/iso_spec_model_ssh.eps',bbox_inches='tight')
plt.savefig('figs/iso_spec_model_ssh.pdf',bbox_inches='tight')
Etot = np.load("Eiso_uv_small_100m.npz")
Ed = np.load("Eiso_uv_small_daily_average_100m.npz")
#GM = np.load("../GM/gm_Ctot.npz")
GM = io.loadmat('/Users/crocha/Dropbox/research/comparisons/GarrettMunkMatlab-master/GM_southern-ocean_npf.mat',squeeze_me=True,struct_as_record=False)
# E0 is the level of energy in the GM specrtrum
# somehow arbitrary, but according to Munk&Garrent
# maybe universal within a factor of 2. Here
# twice the orignal value seems to work best.
E0 = 1.4
kgm, Egm = GM['k'], GM['EGM']
fgm = (kgm>.5e-3)&(kgm<1.)
kgm, Egm = kgm[fgm], E0*Egm[fgm]
kiso, Eiso = Etot['K_w'], Etot['Ew']
kd, Ed = Ed['K_w'], Ed['Ew']
fi = (1./kiso > 5.)
kiso, Eiso = kiso[fi], Eiso[fi]
kd, Ed = kd[fi], Ed[fi]
fi = (1./kiso < 450.)
kiso, Eiso = kiso[fi], Eiso[fi]
kd, Ed = kd[fi], Ed[fi]
Eld = Eiso-Ed
lw,aph=3,.36
fig = plt.figure(facecolor='w', figsize=(12.,10.))
ax1 = fig.add_subplot(111)
ax1.loglog(kiso,Eiso,color=color1,linewidth=lw,
label=r'hourly')
ax1.loglog(kd,Ed,'--',color=color1,linewidth=lw,label=r'daily-averaged')
ax1.loglog(kd,Eld,color=color2,linewidth=lw,label=r'residual')
ax1.loglog(kgm,Egm,'--', color='k',linewidth=2.,alpha=.7)
ax1.text(0.1, .3e-2,u'$GM$')
plt.xlabel('Wavenumber [cpkm]')
plt.ylabel(u'KE spectral density [ m$^{2}$ s$^{-2}$/ cpkm]')
#plt.text(0.7, 4.5, "a", size=35, rotation=0.)
plt.text(1./20., 3., "llc 4320 simulation", size=25, rotation=0.,
ha="center", va="center",
bbox = dict(boxstyle="round",ec='k',fc='w'))
lg = plt.legend(loc=3,title='isotropic spectrum, 1000 m', numpoints=1,ncol=2)
lg.draw_frame(False)
plt.axis((1./1.e3,1./4.,.5/1.e4,5e0))
plt.text(1./5, 2.8, "f", size=35, rotation=0.)
add_second_axis(ax1)
plt.savefig('figs/iso_spec_model_100m_gm',bbox_inches='tight')