A look at the effect of reducing the amount of slip on the tides.
- 5 day runs with all tidal constituents and flux corrected tides. (double corrected for M2).
- control - rn_shlat = 0.5
- partial1 - rn_shlat = 0.1
Note: free slip when rn_shlat = 0 and no slip when rn_shlat = 2.
Purpose: determine if modifying the amount of slip can affect the tidal amplitudes and phases.
%matplotlib inline
import netCDF4 as NC
from salishsea_tools import tidetools
import matplotlib.pyplot as plt
import math
import numpy as np
# grid
bathy, X, Y = tidetools.get_SS_bathy_data()
path = '/data//nsoontie/MEOPAR/SalishSea/results/partial_slip/'
name1 = 'control'
name2 = 'partial1'
location = ("YorkeIsland", "KelseyBay", "BrownBay", "MaudeIslandE", "CampbellRiver", "ChathamPoint", "TwinIslets", "Lund")
Hmm, I'm not sure that we can trust the tidal harmonics output from NEMO since it is only a 5 day run. But I can do a quick plot of the ssh to see if there is any change between runs.
ax = [4240,4380,-2.5,2.5]
plt.figure(figsize=(20,9))
for stn in range(8):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
time1 = fT1.variables["time_counter"][:]/3600. # want hours not seconds
ssh1 = fT1.variables["sossheig"][:,0,0]
fT2 = NC.Dataset(path + name2+ '/' +location[stn]+'.nc','r')
time2 = fT2.variables["time_counter"][:]/3600. # want hours not seconds
ssh2 = fT2.variables["sossheig"][:,0,0]
plt.subplot(3,3,stn+1)
plt.plot (time1, ssh1,label=name1)
plt.plot (time2, ssh2,label=name2)
plt.title(location[stn])
plt.axis(ax)
plt.grid(True)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
plt.figure(figsize=(8,10))
plt.pcolormesh(X,Y,bathy)
for stn in range(8):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
lat = fT1.variables["nav_lat"]
lon = fT1.variables["nav_lon"]
plt.plot(lon[0,0],lat[0,0],'o',label=location[stn])
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
plt.axis([-126.2,-124.5, 49.8,50.6])
[-126.2, -124.5, 49.8, 50.6]
Ok, there seems to be larger elevations in the Johnstone Strait when we have more friction through the lateral BCs, that is the control run has a larger range in SSH at least for Yorke Island and
Very rough estimate of amplitudes: max elevation - minimum elevation
(does it make sense to do this?)
print 'Station / Max-Min(Control) / Max-Min(Partial1)'
for stn in range(8):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
time1 = fT1.variables["time_counter"][:]/3600. # want hours not seconds
ssh1 = fT1.variables["sossheig"][:,0,0]
max1 = np.max(ssh1); min1 = np.min(ssh1)
diff1 = max1 - min1
fT2 = NC.Dataset(path + name2+ '/' +location[stn]+'.nc','r')
time2 = fT2.variables["time_counter"][:]/3600. # want hours not seconds
ssh2 = fT2.variables["sossheig"][:,0,0]
max2 = np.max(ssh2); min2 = np.min(ssh2)
diff2 = max2-min2
print location[stn], ' / ' , diff1, ' / ', diff2
Station / Max-Min(Control) / Max-Min(Partial1) YorkeIsland / 4.23757 / 4.17186 KelseyBay / 4.22347 / 4.16643 BrownBay / 3.63307 / 3.63674 MaudeIslandE / 2.50328 / 2.65081 CampbellRiver / 3.11669 / 3.11681 ChathamPoint / 3.86258 / 3.80611 TwinIslets / 3.74823 / 3.83391 Lund / 3.76178 / 3.84686
There is a reduction in max-min at Yorke Island and Kelsey Bay and Chatham Point. Closer to free slip = smaller amplitudes at these stations, which is what we needed.
Campbell River and Brown Bay are not much affected.
Maude Island E sees the biggest difference using this diagnostic. The difference suggests increased amplitudes at this point.
Surface currents next
U and V are not exactly aligned at the same point. I may consider averaging to a Tgrid point.
This may not be helpful since data is so coarse in time. (6 hour averages)
icoords=[]
jcoords=[]
ax=ax = [4240,4380,0,2]
print 'Currents'
plt.figure(figsize=(15,9))
for stn in range(6):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
lat = fT1.variables["nav_lat"][0,0]
lon = fT1.variables["nav_lon"][0,0]
i,j = tidetools.find_closest_model_point(lon,lat,X,Y,bathy)
icoords.append(i)
jcoords.append(j)
fU1 = NC.Dataset(path + name1+ '/' +'SalishSea_6h_20030421_20030425_grid_U.nc','r')
fV1 = NC.Dataset(path + name1+ '/' +'SalishSea_6h_20030421_20030425_grid_V.nc','r')
fU2 = NC.Dataset(path + name2+ '/' +'SalishSea_6h_20030421_20030425_grid_U.nc','r')
fV2 = NC.Dataset(path + name2+ '/' +'SalishSea_6h_20030421_20030425_grid_V.nc','r')
for stn in range(6):
time = fU1.variables["time_counter"][:]/3600.
u1 = fU1.variables['vozocrtx'][:,0,icoords[stn],jcoords[stn]]
v1 = fV1.variables['vomecrty'][:,0,icoords[stn],jcoords[stn]]
u2 = fU2.variables['vozocrtx'][:,0,icoords[stn],jcoords[stn]]
v2 = fV2.variables['vomecrty'][:,0,icoords[stn],jcoords[stn]]
plt.subplot(2,3,stn+1)
curr1 = np.sqrt(u1**2+v1**2)
plt.plot (time, curr1,label=name1)
curr2 = np.sqrt(u2**2+v2**2)
plt.plot (time, curr2,label=name2)
plt.title(location[stn])
plt.axis(ax)
plt.grid(True)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
plt.figure(figsize=(8,10))
plt.pcolormesh(X,Y,bathy)
for stn in range(6):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
lat = fT1.variables["nav_lat"]
lon = fT1.variables["nav_lon"]
plt.plot(lon[0,0],lat[0,0],'o',label=location[stn])
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
plt.axis([-126.2,-124.8, 49.8,50.6])
Currents
[-126.2, -124.8, 49.8, 50.6]
Yorke Island Currents are zero at Yorke Island. We are probably hitting a masked point on the U/V grids here. Oh well
** What about the rest of the Strait?**
location = ("PedderBay", "PointGrey", "PortTownsend", "HalfmoonBay", "PortRenfrew")
ax = [4240,4380,-2.5,2.5]
plt.figure(figsize=(20,9))
for stn in range(5):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
time1 = fT1.variables["time_counter"][:]/3600. # want hours not seconds
ssh1 = fT1.variables["sossheig"][:,0,0]
fT2 = NC.Dataset(path + name2+ '/' +location[stn]+'.nc','r')
time2 = fT2.variables["time_counter"][:]/3600. # want hours not seconds
ssh2 = fT2.variables["sossheig"][:,0,0]
plt.subplot(2,3,stn+1)
plt.plot (time1, ssh1,label=name1)
plt.plot (time2, ssh2,label=name2)
plt.title(location[stn])
plt.axis(ax)
plt.grid(True)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
plt.figure(figsize=(8,10))
plt.pcolormesh(X,Y,bathy)
for stn in range(5):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
lat = fT1.variables["nav_lat"]
lon = fT1.variables["nav_lon"]
plt.plot(lon[0,0],lat[0,0],'o',label=location[stn])
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
#plt.axis([-126.2,-124.8, 49.8,50.6])
<matplotlib.legend.Legend at 0x6314890>
print 'Station / Max-Min(Control) / Max-Min(Partial1)'
for stn in range(5):
fT1 = NC.Dataset(path + name1+ '/' +location[stn]+'.nc','r')
time1 = fT1.variables["time_counter"][:]/3600. # want hours not seconds
ssh1 = fT1.variables["sossheig"][:,0,0]
max1 = np.max(ssh1); min1 = np.min(ssh1)
diff1 = max1 - min1
fT2 = NC.Dataset(path + name2+ '/' +location[stn]+'.nc','r')
time2 = fT2.variables["time_counter"][:]/3600. # want hours not seconds
ssh2 = fT2.variables["sossheig"][:,0,0]
max2 = np.max(ssh2); min2 = np.min(ssh2)
diff2 = max2-min2
print location[stn], ' / ' , diff1, ' / ', diff2
Station / Max-Min(Control) / Max-Min(Partial1) PedderBay / 2.18453 / 2.17409 PointGrey / 3.49373 / 3.58613 PortTownsend / 2.63806 / 2.64608 HalfmoonBay / 3.61464 / 3.7082 PortRenfrew / 2.59348 / 2.58552
Not much difference. Maybe slightly larger amplitudes with more free slip.
Curious to see what the currents look like up here in the north.
fU1 = NC.Dataset(path + name1+ '/' +'SalishSea_6h_20030421_20030425_grid_U.nc','r')
fV1 = NC.Dataset(path + name1+ '/' +'SalishSea_6h_20030421_20030425_grid_V.nc','r')
fU2 = NC.Dataset(path + name2+ '/' +'SalishSea_6h_20030421_20030425_grid_U.nc','r')
fV2 = NC.Dataset(path + name2+ '/' +'SalishSea_6h_20030421_20030425_grid_V.nc','r')
t = 8
plt.figure(figsize=(20,9))
u1 = fU1.variables['vozocrtx'][t,0,:,:]
v1 = fV1.variables['vomecrty'][t,0,:,:]
u2 = fU2.variables['vozocrtx'][t,0,:,:]
v2 = fV2.variables['vomecrty'][t,0,:,:]
u1=np.ma.masked_values(u1, 0);
u2=np.ma.masked_values(u2, 0);
v1=np.ma.masked_values(v1, 0);
v2=np.ma.masked_values(v2, 0);
ax=[115,125,769,781]
plt.subplot(2,2,1)
plt.pcolormesh(u1,vmin=-1,vmax=1); plt.title('U1'); plt.grid(True)
plt.axis(ax)
plt.subplot(2,2,2)
plt.pcolormesh(v1,vmin=-1,vmax=1); plt.title('V1'); plt.grid(True)
plt.axis(ax)
plt.subplot(2,2,3)
plt.pcolormesh(u2,vmin=-1,vmax=1); plt.title('U2'); plt.grid(True)
plt.axis(ax)
plt.subplot(2,2,4)
plt.pcolormesh(v2,vmin=-1,vmax=1); plt.title('V2'); plt.grid(True)
plt.axis(ax)
[115, 125, 769, 781]
The ugrid looks funny to me. What does T look like?
ax=[115,125,769,781]
plt.pcolormesh(bathy); plt.title('bathy')
plt.axis(ax)
[115, 125, 769, 781]
Summary
I can't actually measure the tidal amplitudes since the run wasn't long enough and I used all the tidal components. But it looks like pushing closer to free slip will decrease amplitudes in the far north and might not affect much of the southern Strait.