Notebook to investigate the strength of fronts associated with the Fraser River Plume
%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import netCDF4 as NC
import numpy as np
from salishsea_tools import viz_tools
grid = NC.Dataset('../../NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
bathy = grid.variables['Bathymetry'][:]
time1 = 0
# Tracer files 1
fT1 = NC.Dataset('/data/sallen/MEOPAR/SalishSea/NorthTides4/5mn_grid_T.nc','r')
sal1 = fT1.variables['vosaline'][time1,:]
depth = fT1.variables['deptht'][:]
lats = fT1.variables['nav_lat'][:]
lons = fT1.variables['nav_lon'][:]
time2 = 0
# Tracer files 2
fT2 = NC.Dataset('/data/sallen/MEOPAR/SalishSea/Sep18_1min/1mn_grid_T.nc','r')
sal2 = fT2.variables['vosaline'][time2,:]
# depth thicknesses
deldepth = np.zeros(40)
deldepth[0] = 0.5*(depth[0]+depth[1])
deldepth[1:39] = 0.5*(depth[1:39]+depth[2:40])-0.5*(depth[0:38]+depth[1:39])
# mask salinity arrays
m = sal1 == 0
sal1_masked = np.ma.array(sal1,mask=m)
m = sal2 == 0
sal2_masked = np.ma.array(sal2,mask=m)
level = 0
vs = (15,25,30.5)
fig, axs = plt.subplots(1, 2, figsize=(12,7), sharey=True)
cmap = plt.get_cmap('spectral')
cmap.set_bad('burlywood')
nicecolours = ('white','blue','black')
# first one
viz_tools.set_aspect(axs[0])
mesh = axs[0].pcolormesh(sal1_masked[level],cmap=cmap)
axs[0].set_title('Oct 18, 2002: nu high')
cbar=fig.colorbar(mesh, ax=axs[0])
CS=axs[0].contour(sal1_masked[level],vs, colors=nicecolours)
cbar.add_lines(CS)
# second one
viz_tools.set_aspect(axs[1])
mesh = axs[1].pcolormesh(sal2_masked[level],cmap=cmap)
axs[1].set_title('Sep 18, 2003: nu low')
cbar=fig.colorbar(mesh, ax=axs[1])
CS=axs[1].contour(sal2_masked[level],vs, colors=nicecolours)
cbar.add_lines(CS)
# a diagonal, approximately the Duke Island Ferry Route
count = 0
ist = 495
ied = 405
di = np.abs(ied-ist+1)
jst = 225
delta = 1.
dl = np.sqrt((di*440)**2+(di*delta*500)**2)/1000./di
points1 =np.zeros(di+1);lats1 = np.zeros(di+1);lons1 = np.zeros(di+1)
points2 = np.zeros(di+1)
for i in range(ist,ied,-1):
j = jst+count*delta
points1[count] = sal1_masked[0,i,j]
lats1[count] = i
lons1[count] = j
points2[count] = sal2_masked[0,i,j]
count = count+1
fig, ax = plt.subplots(1, 2, figsize=(12, 6.5))
ax[0].plot(np.arange(di+1)*dl,points1,'b',label="Oct 2002, High Viscosity")
ax[0].plot(np.arange(di+1)*dl,points2,'m-',label="Sep 2003, Low Viscosity")
ax[0].legend()
viz_tools.set_aspect(ax[1])
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax[1].pcolormesh(bathy, cmap=cmap)
fig.colorbar(mesh)
plt.axis((200, 350, 350, 550))
ax[1].plot(lons1,lats1,'-k')
/home/sallen/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:3847: UserWarning: Warning: converting a masked element to nan. warnings.warn("Warning: converting a masked element to nan.")
[<matplotlib.lines.Line2D at 0x5d216d0>]
# across the domain
i = 405
jst=245
jed=320
dj = np.abs(jed-jst)
dl = 0.44
fig, ax = plt.subplots(1, 2, figsize=(12, 6.5))
ax[0].plot(np.arange(dj)*dl,sal1_masked[0,i,jst:jed],'b',label="Oct 2002, High Viscosity")
ax[0].plot(np.arange(dj)*dl,sal2_masked[0,i,jst:jed],'m-',label="Sep 2003, Low Viscosity")
ax[0].legend()
viz_tools.set_aspect(ax[1])
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax[1].pcolormesh(bathy, cmap=cmap)
fig.colorbar(mesh)
plt.axis((200, 350, 350, 550))
ax[1].plot(np.arange(jst,jed),np.ones(dj)*i,'-k')
[<matplotlib.lines.Line2D at 0x9642c90>]
# along the domain
ist = 380
ied = 480
j = 305
di = np.abs(ied-ist)
dl = 0.5
fig, ax = plt.subplots(1, 2, figsize=(12, 6.5))
ax[0].plot(np.arange(di)*dl,sal1_masked[0,ist:ied,j],'b',label="Oct 2002, High Viscosity")
ax[0].plot(np.arange(di)*dl,sal2_masked[0,ist:ied,j],'m-',label="Sep 2003, Low Viscosity")
ax[0].legend()
viz_tools.set_aspect(ax[1])
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax[1].pcolormesh(bathy, cmap=cmap)
fig.colorbar(mesh)
plt.axis((200, 350, 350, 550))
ax[1].plot(np.ones(di)*j,np.arange(ist,ied),'-k')
[<matplotlib.lines.Line2D at 0xaef6f10>]
# profile
i = 455
j = 280
dd = 15
fig, ax = plt.subplots(1, 2, figsize=(12, 6.5))
ax[0].plot(sal1_masked[:dd,i,j],-depth[:dd],'b+-',label="Oct 2002, High Viscosity")
ax[0].plot(sal2_masked[:dd,i,j],-depth[:dd],'m+-',label="Sep 2003, Low Viscosity")
ax[0].legend(loc="lower left")
viz_tools.set_aspect(ax[1])
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax[1].pcolormesh(bathy, cmap=cmap)
fig.colorbar(mesh)
plt.axis((200, 350, 350, 550))
ax[1].plot(j,i,'k*')
#plt.plot(sal1_masked[:15,455,280],-depth[:15],'bx-',sal2_masked[:15,455,280],-depth[:15],'mx-')
[<matplotlib.lines.Line2D at 0xc4403d0>]
Mark and Rich post the following, (http://www.oceannetworks.ca/fraser-river-plume)
showing an extremely strong front. However, in Halverson and Pawlowicz (2008): http://onlinelibrary.wiley.com/doi/10.1029/2008JC004844/full, they produce a typical profile:
Compared to the latter our fronts are of similar strength.