This notebook is studying the vortices shed by the WCR. Inspired by output from runew-06-spall-const-5d-4-shc

Hypotheses are:

Barotropic instability of near shelfbreak jet - solve the barotropic instability problem for a Gaussian jet.
WCR exports vorticity in bursts
Vortex rollup of sheet exported by WCR



In [1]:

    
%matplotlib inline

from numpy import *
import h5py
import netCDF4
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['savefig.dpi'] = 2 * matplotlib.rcParams['savefig.dpi']
import numexpr as ne
from scipy.io import loadmat

import pyroms









    



scrip.so not found. Remapping function will not be available

First, bring up plots



In [2]:

    
dirname = 'runs/topoeddy/runew-06-spall-const-5d-4-shc/'
fname = dirname +  'ocean_avg.nc.new1'
iname = dirname + '/config/tea_spall_ini.nc'
edname = dirname + 'eddytrack.mat'

# This time index looks good
tind = 20

# fftconvolve is pretty slow...
def avg1(variable, axis):
    import scipy
    from scipy import signal
    cnv = scipy.ones(variable.ndim)
    cnv[axis] = 2
    cnv = scipy.ones(cnv)
    out = signal.fftconvolve(variable, cnv/2, mode='valid')
    return out

# open file
hisfile = h5py.File(fname,'r')
inifile = netCDF4.Dataset(iname,'r')

# get roms grid using pyroms
rgrid = pyroms.grid.get_ROMS_grid('id',hist_file=fname,grid_file=fname)
dz = diff(rgrid.vgrid.z_w[:,:,:], axis=0)

# get variables I need
ubar = hisfile['ubar'][tind,:,:]
vbar = hisfile['vbar'][tind,:,:]
beta = inifile.variables['phys.beta'][:]

# depth integrate - doesn't work because rvor is incorrect in the first place
#rvorbar = sum(rvor * avg1(avg1(dz/rgrid.vgrid.h,1),2), axis=0)

# Calculate rvorbar terms
Uy = diff(ubar,n=1,axis=0)
Vx = diff(vbar,n=1,axis=1)
rvorbar = Vx - Uy;









    



Load cartesian grid from file

Vorticity Export

Plot area-averaged vorticity exported onto the shelf and look at the pattern. I need to calculate pvflux = $ \frac{1}{H} \frac{1}{L_x}\int \int_{-H}^{\zeta} v q \, dz \, dx $ where, q is Ertel PV.



In [3]:

    
isb = 20 ; # must add code to figure this out
hsb = hisfile['h'][26,1]
Lx = amax(rgrid.hgrid.x_rho[:])

# get first eddy shedding
fname2 = dirname + 'ocean_avg.nc'
pvname = dirname + 'ocean_vor.nc'

eddy = loadmat(edname)['eddy']

trange=range(10,65)

zr = rgrid.vgrid.z_r[:]
zr = zr[:,0:isb+5,:]

pvfile = h5py.File(pvname,'r');

dV = avg1(1500 * 1500 * dz,0);
v = hisfile['v'][trange,:,0:isb+5,1:-1]
pv = pvfile['pv'][trange,:,0:isb+5,:]
rv = avg1(pvfile['rv'][trange,:,0:isb+6,:],3)

tpv = pvfile['ocean_time'][trange]

rvflux = sum(sum( rv[:,:,isb,:] * avg1(v[:,:,isb,:]*dz[:,isb,1:-1],1),axis=1),axis=1)/hsb/Lx
pvflux = sum(sum( pv[:,:,isb,:] * avg1(v[:,:,isb,:],1)*avg1(dz[:,isb,1:-1],0),axis=1),axis=1)/hsb/Lx

shelfvol = sum(dV[:,0:isb,:][:])
pvshelf = sum(sum(sum( pv[:,:,0:isb,:] * dV[:,0:isb,1:-1], axis=1), axis=1), axis=1)/shelfvol
rvshelf = sum(sum(sum( rv[:,:,0:isb,:] * dV[:,0:isb,1:-1], axis=1), axis=1), axis=1)/shelfvol



In [73]:

    
plt.subplot(2,1,1)
plt.plot(tpv/86400,pvflux/amax(abs(pvflux[:])),tpv/86400,rvflux/amax(abs(rvflux[:])))
plt.xlabel('Time (days)')
plt.ylabel('Area - averaged fluxes')
plt.legend( ('Ertel PV','Rel Vor.'), loc='lower left')

plt.subplot(2,1,2)
plt.plot(tpv/86400,pvshelf - pvshelf[0], tpv/86400,rvshelf - rvshelf[0])
plt.ylabel('Volume average vorticity on shelf')
plt.xlabel('Time (days)');

Barotropic Instability



In [6]:

    
plt.subplot(2,1,1)
plt.contourf(-Uy[0:100,:],30)
plt.title('U_y')
plt.colorbar()

plt.subplot(2,1,2)
plt.contourf((Vx)[0:100,:])
plt.title('V_x')
plt.colorbar()









    Out[6]:





<matplotlib.colorbar.Colorbar instance at 0x5289050>

Now compare $U_{yy}$ to $\beta$



In [7]:

    
Uyy = diff(ubar,n=2,axis=0)/1500**2
plt.contourf((beta - Uyy)[0:100,:])

plt.colorbar()









    Out[7]:





<matplotlib.colorbar.Colorbar instance at 0x51583b0>



In [337]:

    
plt.plot(Vx[:,100])
plt.hold()
plt.plot(Uy[:,100])
plt.plot( (beta - Uyy[:,100]).transpose()*1e7 )
plt.plot(ubar[:,100]*0.1)
plt.xticks(list(plt.xticks()[0]) + [26,59])
plt.legend( ('V_x','U_y','beta - U_yy*1e7','U*0.1') )
plt.xlabel('Y co-ordinate')









    Out[337]:





<matplotlib.text.Text at 0x497cafd0>

Line plots of $U_y, \beta - U_{yy}, V_x$ . As expected, $V_x \ll U_y$, and so, I can idealize this as a barotropic along-shore jet that is pretty much centered over the slope. Shelfbreak is at 26 and slopebreak at 59

PV gradient $\beta - U_{y}$ crosses zero - Rayleigh's criterion is satisfied. The vorticity maximum is also within the domain, and not at the boundary or inifnity - hence, Fjortoft's criterion is also satisfied.

Barotropic instability is possible! Note that $\beta$ cannot be ignored.

For the theoretical derivation, let's assume that I have a Gaussian velocity profile



In [325]:

    
y = linspace(-pi/2,pi/2,120)

# parameters
U0 = -1
L = 5

# functional forms
U = cos(2*pi*y)**2 #U0 * exp(-(y/L)**2);
Uy = diff(U,n=1)/diff(y)
Uyy = diff(U,n=2)/avg1(diff(y)**2,0)
y1 = avg1(y,0)
y2 = avg1(y1,0)

# plot
plt.plot(y,U, 'c', y1,Uy,'g', y2,-Uyy,'r')
plt.legend( ('U','Uy','-Uyy') )









    Out[325]:





<matplotlib.legend.Legend at 0x41ddb750>

I need an analytical expression for the growth rate now. Then code it. Pick out appropriate values for the Gaussian assumption and then see if the estimates match.



In [ ]: