Make tidal boundary conditions for the Northern boundary in Johnstone Strait. Use data from Webtide and from Thomson & Huggett (1980) to make netcdf input files.
Edited by Nancy in March, 2014. Modified measured grid rotation ($129^\circ$ instead of $-51^\circ$) and error in S2 elevation calculation.
Edited by Susan in April, 2014. Modified grid rotation to rotate versus across grid rather than up grid. Corrects phases of currents.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import webtide_tools
In [2]:
TSin = {'North': {'M2': {'Amp': 1.2468500142849095, 'Pha': 5.8718032013900494}, 'Q1': {'Amp': 0.88018411744300051, 'Pha': 3.3676909436466511}, 'O1': {'Amp': 0.92187588787164776, 'Pha': -1.8419422569093058}, 'S2': {'Amp': 0.98646859951016619, 'Pha': -10.363961743466264}, 'K1': {'Amp': 1.056369887585026, 'Pha': -7.2574882761515624}, 'N2': {'Amp': 0.95930333678587432, 'Pha': -6.8740655239509989}}, 'West': {'M2': {'Amp': 0.8840061749307605, 'Pha': -7.8587198373167935}, 'Q1': {'Amp': 0.95936152535045738, 'Pha': 1.1445537519586795}, 'O1': {'Amp': 1.058133345863022, 'Pha': 0.92064288566727737}, 'S2': {'Amp': 1.0232932221404019, 'Pha': 0.7966194737248522}, 'K1': {'Amp': 1.0168859739167362, 'Pha': -3.8356999030330829}, 'N2': {'Amp': 1.0027750731988296, 'Pha': 2.5899670236940464}}}
In [3]:
# These are the corrections to be used.
tidesetname = 'N36_AF'
M2_amp_finecorrection = TSin['North']['M2']['Amp']
M2_pha_finecorrection = TSin['North']['M2']['Pha']
K1_amp_finecorrection = TSin['North']['K1']['Amp']
K1_pha_finecorrection = TSin['North']['K1']['Pha']
In [4]:
#from Table 2, Thomson & Huggett (1980)
#Tidal height harmonics
#lats and lons found by estimating locations shown in Fig 4 with Google Maps
TGS_Jun = {'O1amp':0.3249,'O1pha_PST':130.3,
'K1amp':0.5491,'K1pha_PST':138.9,
'M2amp':1.2196,'M2pha_PST':37.3,
'S2amp':0.3966,'S2pha_PST':56.6,
'Lat':50.470945,'Lon':-126.23826}
TGN_Jun = {'O1amp':0.3250,'O1pha_PST':129.4,
'K1amp':0.5504,'K1pha_PST':138.5,
'M2amp':1.2086,'M2pha_PST':36.8,
'S2amp':0.3931,'S2pha_PST':56.6,
'Lat':50.498579,'Lon':-126.213026}
TGW_Jan = {'O1amp':0.3161,'O1pha_PST':127.5,
'K1amp':0.5466,'K1pha_PST':138.7,
'M2amp':1.2453,'M2pha_PST':28.8,
'S2amp':0.4504,'S2pha_PST':51.2,
'Lat':50.574707,'Lon':-126.790431}
TGS_Feb = {'O1amp':0.3198,'O1pha_PST':126.8,
'K1amp':0.5696,'K1pha_PST':139.4,
'M2amp':1.2218,'M2pha_PST':37.5,
'S2amp':0.4474,'S2pha_PST':59.1,
'Lat':50.470945,'Lon':-126.23826}
TGE_Feb = {'O1amp':0.3376,'O1pha_PST':127.5,
'K1amp':0.5881,'K1pha_PST':140.8,
'M2amp':1.1690,'M2pha_PST':41.1,
'S2amp':0.4368,'S2pha_PST':62.1,
'Lat':50.442638,'Lon':-125.97064}
In [5]:
#need to convert phases from relative to 120 degrees W to GMT
#for this conversion, we need to use the frequency of each component (in hours)
freqs = {'M2':12.42, 'S2':12.00, 'O1':25.82, 'K1':23.93}
In [6]:
#For details, see Manual for Tide Heights Analysis and Prediction by M.G.G. Foreman,
#Pacific Marine Science Report 77-10, IOS, 1977 (Revised 2004),
#Section 2.3.1 Astronomical argument and Greenwich phase lag
# available here: ftp://canuck.seos.uvic.ca/docs/MFTides/heights.pdf
# or here: http://www.pac.dfo-mpo.gc.ca/science/oceans/tidal-marees/index-eng.html
#Now, 120 degrees W (Pacific Standard Time) is 8 hours different to GMT
#So to convert M2 phase from PST to GMT:
# GMT_pha = PST_pha - (-8hours)*(1/12.42hours)*360degrees
#TGS_Jun
TGS_Jun['O1pha_GMT'] = TGS_Jun['O1pha_PST']-(-8.0)*(1.0/freqs['O1'])*360
TGS_Jun['K1pha_GMT'] = TGS_Jun['K1pha_PST']-(-8.0)*(1.0/freqs['K1'])*360
TGS_Jun['M2pha_GMT'] = TGS_Jun['M2pha_PST']-(-8.0)*(1.0/freqs['M2'])*360
TGS_Jun['S2pha_GMT'] = TGS_Jun['S2pha_PST']-(-8.0)*(1.0/freqs['S2'])*360
#TGN_Jun
TGN_Jun['O1pha_GMT'] = TGN_Jun['O1pha_PST']-(-8.0)*(1.0/freqs['O1'])*360
TGN_Jun['K1pha_GMT'] = TGN_Jun['K1pha_PST']-(-8.0)*(1.0/freqs['K1'])*360
TGN_Jun['M2pha_GMT'] = TGN_Jun['M2pha_PST']-(-8.0)*(1.0/freqs['M2'])*360
TGN_Jun['S2pha_GMT'] = TGN_Jun['S2pha_PST']-(-8.0)*(1.0/freqs['S2'])*360
#TGW_Jan
TGW_Jan['O1pha_GMT'] = TGW_Jan['O1pha_PST']-(-8.0)*(1.0/freqs['O1'])*360
TGW_Jan['K1pha_GMT'] = TGW_Jan['K1pha_PST']-(-8.0)*(1.0/freqs['K1'])*360
TGW_Jan['M2pha_GMT'] = TGW_Jan['M2pha_PST']-(-8.0)*(1.0/freqs['M2'])*360
TGW_Jan['S2pha_GMT'] = TGW_Jan['S2pha_PST']-(-8.0)*(1.0/freqs['S2'])*360
#TGS_Feb
TGS_Feb['O1pha_GMT'] = TGS_Feb['O1pha_PST']-(-8.0)*(1.0/freqs['O1'])*360
TGS_Feb['K1pha_GMT'] = TGS_Feb['K1pha_PST']-(-8.0)*(1.0/freqs['K1'])*360
TGS_Feb['M2pha_GMT'] = TGS_Feb['M2pha_PST']-(-8.0)*(1.0/freqs['M2'])*360
TGS_Feb['S2pha_GMT'] = TGS_Feb['S2pha_PST']-(-8.0)*(1.0/freqs['S2'])*360
#TGE_Feb
TGE_Feb['O1pha_GMT'] = TGE_Feb['O1pha_PST']-(-8.0)*(1.0/freqs['O1'])*360
TGE_Feb['K1pha_GMT'] = TGE_Feb['K1pha_PST']-(-8.0)*(1.0/freqs['K1'])*360
TGE_Feb['M2pha_GMT'] = TGE_Feb['M2pha_PST']-(-8.0)*(1.0/freqs['M2'])*360
TGE_Feb['S2pha_GMT'] = TGE_Feb['S2pha_PST']-(-8.0)*(1.0/freqs['S2'])*360
In [7]:
#Compare the harmonics to Webtide harmonics at 50.556435 deg N, -126.7717 deg W
import pandas
#webtide = pandas.read_csv('/ocean/klesouef/meopar/tools/I_ForcingFiles/Tides/webtide_point_johnstone_strait.csv',
#index_col=0,skiprows=2)
webtide = pandas.read_csv('webtide_point_johnstone_strait.csv',
index_col=0,skiprows=2)
webtide = webtide.rename(columns={'Amplitude (m)':'amp','Phase (deg GMT)':'pha'})
print (webtide)
In [8]:
#Plot them to check how it looks
#the channel is mostly aligned east-west, so plot along the channel
plt.figure(figsize=[9,5])
pha = 'M2pha_GMT','O1pha_GMT','K1pha_GMT','S2pha_GMT'
cols = 'r','b','m','k','g'
marks = '>','^','*','d'
for t in range(0,len(pha)):
plt.plot([TGW_Jan['Lon'],TGS_Feb['Lon'],TGS_Jun['Lon'],TGN_Jun['Lon'],TGE_Feb['Lon']],
[TGW_Jan[pha[t]],TGS_Feb[pha[t]],TGS_Jun[pha[t]],TGN_Jun[pha[t]],TGE_Feb[pha[t]]],
cols[t],marker=marks[t],markersize=8,linestyle='none')
pha = 'M2','O1','K1','S2'
nums = [6,0,3,7]
marknum = 0
for t in nums:
plt.plot(webtide['Longitude'][t],webtide.pha[t],color='green',marker=marks[marknum],markersize=8,linestyle='none')
marknum = marknum+1
plt.xlim((-127, -125.5))
plt.xlabel('Longitude (deg W)')
plt.ylabel('Phase (deg GMT)')
plt.legend(('M2pha_GMT','O1pha_GMT','K1pha_GMT','S2pha_GMT'),loc=4)
plt.title('Phase lag for measured and Webtide points')
Out[8]:
In [9]:
#Plot them to check how it looks
#the channel is mostly aligned east-west, so plot along the channel
plt.figure(figsize=[9,5])
amp = 'M2amp','O1amp','K1amp','S2amp'
cols = 'r','b','m','k','g'
marks = '>','^','*','d'
for t in range(0,len(pha)):
plt.plot([TGW_Jan['Lon'],TGS_Feb['Lon'],TGS_Jun['Lon'],TGN_Jun['Lon'],TGE_Feb['Lon']],
[TGW_Jan[amp[t]],TGS_Feb[amp[t]],TGS_Jun[amp[t]],TGN_Jun[amp[t]],TGE_Feb[amp[t]]],
cols[t],marker=marks[t],markersize=8,linestyle='none')
amp = 'M2','O1','K1','S2'
nums = [6,0,3,7]
marknum = 0
for t in nums:
plt.plot(webtide['Longitude'][t],webtide.amp[t],color='green',marker=marks[marknum],markersize=8,linestyle='none')
marknum = marknum+1
plt.xlim((-127, -125.5))
plt.xlabel('Longitude (deg W)')
plt.ylabel('Amplitude (cm)')
plt.legend(('M2amp','O1amp','K1amp','S2amp'),loc=4)
plt.title('Amplitude for measured and Webtide points')
Out[9]:
In [10]:
#get the measured tidal currents in Johnstone Strait
#table 1 of Thomson & Huggett (1980)
#thomson = pandas.read_csv('/ocean/klesouef/meopar/tools/I_ForcingFiles/Tides/thomson_obs_1980_johnstone_strait.csv')
thomson = pandas.read_csv('thomson_obs_1980_johnstone_strait.csv')
thomson = thomson.rename(columns={'K1 amp (cm/s)':'k1amp','K1 phase (deg relative to 120deg W)':'k1pha_PST',
'M2 amp (cm/s)':'m2amp','M2 phase (deg relative to 120deg W)':'m2pha_PST',})
thomson.k1amp[thomson.k1amp>100] = float('NaN')
thomson.k1pha_PST[thomson.k1pha_PST>1000] = float('NaN')
print (thomson)
In [11]:
#convert the GMT phase to PST phase and rename it to avoid confusion
thomson.k1pha_PST = thomson.k1pha_PST-(-8.0)*(1.0/freqs['K1'])*360
thomson.m2pha_PST = thomson.m2pha_PST-(-8.0)*(1.0/freqs['M2'])*360
thomson = thomson.rename(columns={'k1pha_PST':'k1pha_GMT','m2pha_PST':'m2pha_GMT'})
print (thomson)
Plot the location of all measurements on our grid
In [12]:
import netCDF4 as NC
#grid = NC.Dataset('/ocean/klesouef/meopar/nemo-forcing/grid/bathy_meter_SalishSea2.nc','r')
#grid = NC.Dataset('../../../nemo-forcing/grid/bathy_meter_SalishSea2.nc','r')
grid = NC.Dataset('../../../NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
bathy = grid.variables['Bathymetry']
lat = grid.variables['nav_lat']
lon = grid.variables['nav_lon']
plt.figure(figsize=(12,9))
plt.contourf(lon[:],lat[:],bathy[:],range(4,420))
plt.xlim((-126.8,-125.4))
plt.ylim((50.2,50.6))
#I need the lat/long of each of the current measurement points
#From Google Maps, I estimate the following locations:
#Station 1
plt.plot(-126.133747,50.482197,'m*',markersize=15)
#Station 2
plt.plot(-126.135936,50.478701,'m*',markersize=15)
#Station 3
plt.plot(-126.136451,50.472256,'m*',markersize=15)
#Station 4
plt.plot(-126.137567,50.464225,'m*',markersize=15)
#Station 5
plt.plot(-126.138554,50.458815,'m*',markersize=15)
#Webtide point
plt.plot(-126.7717,50.556435,'g*',markersize=15)
#current measurement points
plt.plot([TGW_Jan['Lon'],TGS_Feb['Lon'],TGS_Jun['Lon'],TGN_Jun['Lon'],TGE_Feb['Lon']],
[TGW_Jan['Lat'],TGS_Feb['Lat'],TGS_Jun['Lat'],TGN_Jun['Lat'],TGE_Feb['Lat']],'r*',markersize=15)
#also mark the edge points
import numpy as np
j=np.array([37,48])
plt.plot(lon[897,j],lat[897,j],'og')
plt.title('Location of measurements around our grid')
#so the current measurements are pretty much exactly at our boundary, YAY!
#and the TGN and TGS stations are pretty close for the water elevation
Out[12]:
Now, we need to get depth integrated currents for putting into NEMO. Assume each current meter represents the water from halfway between it and the current meter above it, down to halfway between it and the current meter below it
In [13]:
#for K1
K1_amp_depth_int = np.zeros((5,1))
K1_pha_depth_int = np.zeros((5,1))
#Station 1: only one point, therefore this represents all depths
K1_amp_depth_int[0,0] = thomson.k1amp[thomson.Station == 1]
K1_pha_depth_int[0,0] = thomson.k1pha_GMT[thomson.Station == 1]
#Station 2: depth integrate amplitude and phase
K1_amp_depth_int[1,0] = np.sum([np.mean(thomson.k1amp[np.logical_and(thomson.Station == 2,thomson.depth == 15)])*45,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 2,thomson.depth == 75)])*67.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 2,thomson.depth == 150)])*75,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 2,thomson.depth == 225)])*72.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 2,thomson.depth == 295)])*40])/300.0
K1_pha_depth_int[1,0]= np.sum([np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 15)])*45,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 75)])*67.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 150)])*75,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 225)])*72.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 295)])*40])/300.0
#Station 3: depth integrate amplitude and phase
K1_amp_depth_int[2,0] = np.sum([np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 15)])*17.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 20)])*4,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 23)])*27.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 75)])*63.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 150)])*62.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 200)])*37.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 225)])*25,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 250)])*25,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 3,thomson.depth == 275)])*37.5])/300.0
K1_pha_depth_int[2,0] = np.sum([np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 15)])*17.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 20)])*4,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 23)])*27.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 75)])*63.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 150)])*62.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 200)])*37.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 225)])*25,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 250)])*25,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 275)])*37.5])/300.0
#Station 4: depth integrate amplitude and phase
K1_amp_depth_int[3,0] = np.sum([np.mean(thomson.k1amp[np.logical_and(thomson.Station == 4,thomson.depth == 75)])*112.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 4,thomson.depth == 150)])*87.5,
np.mean(thomson.k1amp[np.logical_and(thomson.Station == 4,thomson.depth == 250)])*100])/300.0
K1_pha_depth_int[3,0] = np.sum([np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 4,thomson.depth == 75)])*112.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 4,thomson.depth == 150)])*87.5,
np.mean(thomson.k1pha_GMT[np.logical_and(thomson.Station == 4,thomson.depth == 250)])*100])/300.0
#Station 5: only one point, therefore this represents all depths
K1_amp_depth_int[4,0] = thomson.k1amp[thomson.Station == 5]
K1_pha_depth_int[4,0] = thomson.k1pha_GMT[thomson.Station == 5]
#for M2
M2_amp_depth_int = np.zeros((5,1))
M2_pha_depth_int = np.zeros((5,1))
#Station 1: only one point, therefore this represents all depths
M2_amp_depth_int[0,0] = thomson.m2amp[thomson.Station == 1]
M2_pha_depth_int[0,0] = thomson.m2pha_GMT[thomson.Station == 1]
#Station 2: depth integrate amplitude and phase
M2_amp_depth_int[1,0] = np.sum([np.mean(thomson.m2amp[np.logical_and(thomson.Station == 2,thomson.depth == 15)])*45,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 2,thomson.depth == 75)])*67.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 2,thomson.depth == 150)])*75,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 2,thomson.depth == 225)])*72.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 2,thomson.depth == 295)])*40])/300.0
M2_pha_depth_int[1,0]= np.sum([np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 15)])*45,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 75)])*67.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 150)])*75,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 225)])*72.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 2,thomson.depth == 295)])*40])/300.0
#Station 3: depth integrate amplitude and phase
M2_amp_depth_int[2,0] = np.sum([np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 15)])*17.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 20)])*4,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 23)])*27.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 75)])*63.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 150)])*62.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 200)])*37.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 225)])*25,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 250)])*25,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 3,thomson.depth == 275)])*37.5])/300.0
M2_pha_depth_int[2,0] = np.sum([np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 15)])*17.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 20)])*4,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 23)])*27.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 75)])*63.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 150)])*62.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 200)])*37.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 225)])*25,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 250)])*25,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 3,thomson.depth == 275)])*37.5])/300.0
#Station 4: depth integrate amplitude and phase
M2_amp_depth_int[3,0] = np.sum([np.mean(thomson.m2amp[np.logical_and(thomson.Station == 4,thomson.depth == 75)])*112.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 4,thomson.depth == 150)])*87.5,
np.mean(thomson.m2amp[np.logical_and(thomson.Station == 4,thomson.depth == 250)])*100])/300.0
M2_pha_depth_int[3,0] = np.sum([np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 4,thomson.depth == 75)])*112.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 4,thomson.depth == 150)])*87.5,
np.mean(thomson.m2pha_GMT[np.logical_and(thomson.Station == 4,thomson.depth == 250)])*100])/300.0
#Station 5: only one point, therefore this represents all depths
M2_amp_depth_int[4,0] = thomson.m2amp[thomson.Station == 5]
M2_pha_depth_int[4,0] = thomson.m2pha_GMT[thomson.Station == 5]
#measurements were in cm/s, convert to m/s
K1_amp_depth_int = K1_amp_depth_int/100.0
M2_amp_depth_int = M2_amp_depth_int/100.0
plt.figure()
plt.subplot(2,1,1)
plt.plot(range(0,5),M2_amp_depth_int,'.-',range(0,5),K1_amp_depth_int,'.-')
plt.legend(('M2','K1'))
plt.title('Initial depth interpolated currents')
plt.subplot(2,1,2)
plt.plot(range(0,5),M2_pha_depth_int,'.-',range(0,5),K1_pha_depth_int,'.-')
#after plotting them, the phase for K1 at Station 1 looks too high, so just take the value for Station 2
K1_pha_depth_int[0,0] = K1_pha_depth_int[1,0]
plt.figure()
plt.subplot(2,1,1)
plt.plot(range(0,5),M2_amp_depth_int,'.-',range(0,5),K1_amp_depth_int,'.-')
plt.legend(('M2','K1'))
plt.title('Depth interpolated currents after adjusting one point')
plt.subplot(2,1,2)
plt.plot(range(0,5),M2_pha_depth_int,'.-',range(0,5),K1_pha_depth_int,'.-')
Out[13]:
Plot depths along the boundary
In [14]:
k = range(0,60)
plt.plot(k,bathy[897,k])
plt.gca().invert_yaxis()
plt.xlabel('grid point')
plt.ylabel('depth [m]')
plt.title('Depths at northern boundary')
plt.xlim((36,49))
#so the boundary is from j=37 to j=48 at i=897
Out[14]:
In [15]:
# compare the size of the flux through the Thomson section and the size of the flux through
# our section, remembering that our flux is at 39 degrees to the boundary.
# Start with area
print ('Estimate of Thomson x-sectional area in mkm', 0.55*0.5*300 + (2.4-0.55)*300 + 0.37*300*0.5)
dx = 0.4155 # our grid size
dxc = dx*np.cos(39.*np.pi/180.) # corrected for slant of boundary
print ('Our effective x-sectional area in mkm', np.sum(bathy[897,37:48])*dxc)
# Thomson bathymetry
thomx = (0, 0.18, 0.55, 0.8, 1.4, 1.95, 2.4, 2.46, 2.77)/np.cos(0)+0.5
thomy = (0, 100, 300, 300, 300, 300, 300, 150, 0)
plt.plot (thomx,thomy,'x-')
# Model bathymetry, corrected for slant
krange = np.arange(0,49-36+1)*dxc
depth = np.zeros(49-36+1)
depth[0:49-36+1] = bathy[897,36:50]
plt.plot (krange,depth,'*-')
Out[15]:
In [16]:
# Now estimate flux (M2)
obs_flux = 0.
model_flux = 0.
# Station 5
obs_flux = obs_flux + 0.55 * M2_amp_depth_int[4] * 0.5*300
model_flux = model_flux + sum(depth[0:3] * dx * M2_amp_depth_int[4]*np.cos(39.*np.pi/180.))
# Station 4
obs_flux = obs_flux + (1.1-0.55) * M2_amp_depth_int[3] * 300
model_flux = model_flux + sum(depth[3:5] * dx * M2_amp_depth_int[3]*np.cos(39.*np.pi/180.))
# Station 3
obs_flux = obs_flux + (1.7-1.1) * M2_amp_depth_int[2] * 300
model_flux = model_flux + sum(depth[5:7] * dx * M2_amp_depth_int[2]*np.cos(39.*np.pi/180.))
print (M2_amp_depth_int[2])
# Station 2
obs_flux = obs_flux + (2.4-1.7) * M2_amp_depth_int[1] * 300
model_flux = model_flux + sum(depth[7:9] * dx * M2_amp_depth_int[1]*np.cos(39.*np.pi/180.))
# Station 1
obs_flux = obs_flux + (2.77-2.4) * M2_amp_depth_int[0] * 300 * 0.5
model_flux = model_flux + sum(depth[9:13] * dx * M2_amp_depth_int[0]*np.cos(39.*np.pi/180.))
print (obs_flux, model_flux, model_flux/obs_flux, 856/693.)
print (obs_flux/693., model_flux/856.)
M2_corr = obs_flux/model_flux
M2_corr = M2_corr*M2_corr #reduce by an extra factor for a better match with observations
print (M2_corr)
In [17]:
# And flux K1
obs_flux = 0.
model_flux = 0.
# Station 5
obs_flux = obs_flux + 0.55 * K1_amp_depth_int[4] * 0.5*300
model_flux = model_flux + sum(depth[0:3] * dx * K1_amp_depth_int[4]*np.cos(39.*np.pi/180.))
# Station 4
obs_flux = obs_flux + (1.1-0.55) * K1_amp_depth_int[3] * 300
model_flux = model_flux + sum(depth[3:5] * dx * K1_amp_depth_int[3]*np.cos(39.*np.pi/180.))
# Station 3
obs_flux = obs_flux + (1.7-1.1) * K1_amp_depth_int[2] * 300
model_flux = model_flux + sum(depth[5:7] * dx * K1_amp_depth_int[2]*np.cos(39.*np.pi/180.))
# Station 2
obs_flux = obs_flux + (2.4-1.7) * K1_amp_depth_int[1] * 300
model_flux = model_flux + sum(depth[7:9] * dx * K1_amp_depth_int[1]*np.cos(39.*np.pi/180.))
# Station 1
obs_flux = obs_flux + (2.77-2.4) * K1_amp_depth_int[0] * 300 * 0.5
model_flux = model_flux + sum(depth[9:13] * dx * K1_amp_depth_int[0]*np.cos(39.*np.pi/180.))
print (obs_flux, model_flux, model_flux/obs_flux, 856/693.)
K1_corr = obs_flux/model_flux
print (K1_corr)
So... we have current measurements close to the boundary, but not perfectly aligned with the boundary. But without any other well resolved data west of the NEMO boundary, we can't really interpolate to the boundary location. Plus, there isn't much of a change in tidal amplitudes or phases across the boundary. So just assume that the current measurements represent the currents at the boundary without any interpolation or change.
This is from Susan's notebook (JohnstoneStraitBoundary) where she looked at the boundary
In [18]:
#start and end points
starti = 32
endi = 62
lengthi = endi-starti
So let's assume for currents (Kates original assumptions and Susan's suggestions)
Station 1 represents cells 47 - 62 (16 cells) -- SEA 45-49, 50-62 (18 cells)
Station 2 represents cells 44 - 46 (3 cells) -- SEA 43-44 (2 cells)
Station 3 represents cells 41 - 43 (3 cells) --- SEA 41-42 (2 cells)
Station 4 represents cells 39 - 40 (2 cells) --- SEA 39-40 (2 cells)
Station 5 represents cells 32 - 38 (7 cells) --- SEA 32-35, 36-38 (7 cells)
Put all this into some new variables that are the same length as the boundary
This is as far as I have gotten with the corrections. The position of the fluxes, as given by SEA above needs to be corrected. The impact of the correction from Stn 3 to Stn 4 needs to be propagated. The currents need to be downgraded by the amount estimated above.
Nancy's continuation¶
Plan:
- Correct the currents and check.
- Correct grid positioning.
- Propogate S3-S4 changes through.
Current corrections:
We adjust the currents by multiplying by the correction factors give above:
$U_{mod} = corr*U_{obs}$
In [19]:
print (M2_amp_depth_int)
K1_amp_depth_int_corr = K1_corr*K1_amp_depth_int[:];
M2_amp_depth_int_corr = M2_corr*M2_amp_depth_int[:];
print (M2_amp_depth_int_corr)
In [20]:
#Restimate Fluxes
# Now estimate flux (M2)
obs_flux = 0.
model_flux = 0.
# Station 5
obs_flux = obs_flux + 0.55 * M2_amp_depth_int[4] * 0.5*300
model_flux = model_flux + sum(depth[0:3] * dx * M2_amp_depth_int_corr[4]*np.cos(39.*np.pi/180.))
# Station 4
obs_flux = obs_flux + (1.1-0.55) * M2_amp_depth_int[3] * 300
model_flux = model_flux + sum(depth[3:5] * dx * M2_amp_depth_int_corr[3]*np.cos(39.*np.pi/180.))
# Station 3
obs_flux = obs_flux + (1.7-1.1) * M2_amp_depth_int[2] * 300
model_flux = model_flux + sum(depth[5:7] * dx * M2_amp_depth_int_corr[2]*np.cos(39.*np.pi/180.))
# Station 2
obs_flux = obs_flux + (2.4-1.7) * M2_amp_depth_int[1] * 300
model_flux = model_flux + sum(depth[7:9] * dx * M2_amp_depth_int_corr[1]*np.cos(39.*np.pi/180.))
# Station 1
obs_flux = obs_flux + (2.77-2.4) * M2_amp_depth_int[0] * 300 * 0.5
model_flux = model_flux + sum(depth[9:13] * dx * M2_amp_depth_int_corr[0]*np.cos(39.*np.pi/180.))
print (obs_flux, model_flux, model_flux/obs_flux, 856/693.)
print (obs_flux/693., model_flux/856.)
In [21]:
# And flux K1
obs_flux = 0.
model_flux = 0.
# Station 5
obs_flux = obs_flux + 0.55 * K1_amp_depth_int[4] * 0.5*300
model_flux = model_flux + sum(depth[0:3] * dx * K1_amp_depth_int_corr[4]*np.cos(39.*np.pi/180.))
# Station 4
obs_flux = obs_flux + (1.1-0.55) * K1_amp_depth_int[3] * 300
model_flux = model_flux + sum(depth[3:5] * dx * K1_amp_depth_int_corr[3]*np.cos(39.*np.pi/180.))
# Station 3
obs_flux = obs_flux + (1.7-1.1) * K1_amp_depth_int[2] * 300
model_flux = model_flux + sum(depth[5:7] * dx * K1_amp_depth_int_corr[2]*np.cos(39.*np.pi/180.))
# Station 2
obs_flux = obs_flux + (2.4-1.7) * K1_amp_depth_int[1] * 300
model_flux = model_flux + sum(depth[7:9] * dx * K1_amp_depth_int_corr[1]*np.cos(39.*np.pi/180.))
# Station 1
obs_flux = obs_flux + (2.77-2.4) * K1_amp_depth_int[0] * 300 * 0.5
model_flux = model_flux + sum(depth[9:13] * dx * K1_amp_depth_int_corr[0]*np.cos(39.*np.pi/180.))
print (obs_flux, model_flux, model_flux/obs_flux, 856/693.)
The model and observed fluxes match now. Note that there are likely some integration errors in calculating the fluxes but we probably don't have to worry because these are estimates anways.
Now, correct the positioning based on Susn's corrections. Make sure use the new array for corrected depth integrated currents.
Fine Tuning¶
So above is correct theoretically, but now given the model results, lets do some fine tuning. Based on the run tide_M2phase, the average Kelsey Bay/Yorke Island M2 amplitudes are 3.4% high and the phase is 8.2 degrees low. Add these corrections: here to the velocities and below to the elevations.
In [22]:
# Fine tuning corrections at the top of the notebook
M2_amp_depth_int_corr = M2_amp_depth_int_corr * M2_amp_finecorrection
M2_pha_depth_int = M2_pha_depth_int + M2_pha_finecorrection
K1_amp_depth_int_corr = K1_amp_depth_int_corr * K1_amp_finecorrection
K1_pha_depth_int = K1_pha_depth_int + K1_pha_finecorrection
In [23]:
#K1 current amplitude
amp_K1_bound = np.zeros((lengthi,1))
amp_K1_bound[0:7,0] = K1_amp_depth_int_corr[5-1,0]
amp_K1_bound[7:9,0] = K1_amp_depth_int_corr[4-1,0]
amp_K1_bound[9:11,0] = K1_amp_depth_int_corr[3-1,0]
amp_K1_bound[11:13,0] = K1_amp_depth_int_corr[2-1,0]
amp_K1_bound[13:,0] = K1_amp_depth_int_corr[1-1,0]
#K1 current phase
pha_K1_bound = np.zeros((lengthi,1))
pha_K1_bound[0:7,0] = K1_pha_depth_int[5-1,0]
pha_K1_bound[7:9,0] = K1_pha_depth_int[4-1,0]
pha_K1_bound[9:11,0] = K1_pha_depth_int[3-1,0]
pha_K1_bound[11:13,0] = K1_pha_depth_int[2-1,0]
pha_K1_bound[13:,0] = K1_pha_depth_int[1-1,0]
#M2 current amplitude
amp_M2_bound = np.zeros((lengthi,1))
amp_M2_bound[0:7,0] = M2_amp_depth_int_corr[5-1,0]
amp_M2_bound[7:9,0] = M2_amp_depth_int_corr[4-1,0]
amp_M2_bound[9:11,0] = M2_amp_depth_int_corr[3-1,0]
amp_M2_bound[11:13,0] = M2_amp_depth_int_corr[2-1,0]
amp_M2_bound[13:,0] = M2_amp_depth_int_corr[1-1,0]
#M2 current phase
pha_M2_bound = np.zeros((lengthi,1))
pha_M2_bound[0:7,0] = M2_pha_depth_int[5-1,0]
pha_M2_bound[7:9,0] = M2_pha_depth_int[4-1,0]
pha_M2_bound[9:11,0] = M2_pha_depth_int[3-1,0]
pha_M2_bound[11:13,0] = M2_pha_depth_int[2-1,0]
pha_M2_bound[13:,0] = M2_pha_depth_int[1-1,0]
Now, as noted in the text of Thomson & Huggett (1980), the current amplitudes are in the along-strait direction. Across-strait currents are assumed to be negligible. From the direction of the mean currents, the average along-strait direction is approximately 170 degrees/350 degrees, measured counter-clockwise from east.
So, we must convert the along-strait amplitude and phase into u and v amplitude and phase relative for NEMO. This means that we must rotate them because NEMO wants velocities relative to its grid orientation (not North South).
Following the approach used to rotate the Webtide tidal constituents, we can do the following:
- assume the postive along strait direction is approximately 10 degrees measured clockwise from East
- 'up' in NEMO is rotated 29 degrees counter-clockwise of North
- 'across' in NEMO is rotated 29 degrees counter-clockwise of East
- the angle of rotation between along-strait direction and 'across' in NEMO is therefore $\theta = 29^\circ +10^\circ=39^\circ$
- neglect across-strait velocity (i.e. $v_a = 0$, $\phi_v = 0$)
- $u_a$ and $\phi_u$ are the depth averaged amplitude and phase at each station (calculated above)
following the approach described in the docs (docs/tides/tides_data_acquisition):
so for U:
$Z1 = u_a\cos\theta\cos\phi_u$
$Z2 = u_a\cos\theta\sin\phi_u$
and for V:
$Z1 = -u_a\sin\theta\cos\phi_u$
$Z2 = -u_a\sin\theta\sin\phi_u$
In [24]:
theta = 39
#for K1 currents (these were corrected above)
UZ1_K1 = amp_K1_bound*np.cos(np.radians(theta))*np.cos(np.radians(pha_K1_bound))
UZ2_K1 = amp_K1_bound*np.cos(np.radians(theta))*np.sin(np.radians(pha_K1_bound))
VZ1_K1 = -amp_K1_bound*np.sin(np.radians(theta))*np.cos(np.radians(pha_K1_bound))
VZ2_K1 = -amp_K1_bound*np.sin(np.radians(theta))*np.sin(np.radians(pha_K1_bound))
#for M2 currents
UZ1_M2 = amp_M2_bound*np.cos(np.radians(theta))*np.cos(np.radians(pha_M2_bound))
UZ2_M2 = amp_M2_bound*np.cos(np.radians(theta))*np.sin(np.radians(pha_M2_bound))
VZ1_M2 = -amp_M2_bound*np.sin(np.radians(theta))*np.cos(np.radians(pha_M2_bound))
VZ2_M2 = -amp_M2_bound*np.sin(np.radians(theta))*np.sin(np.radians(pha_M2_bound))
print (pha_K1_bound[0])
And let's assume for surface elevation that the whole boundary is represented by the average between TGN and TGS.
Then convert the amplitude and phase into cosine and sine format that NEMO likes.
In [25]:
#K1 surface elevation
elev_amp_K1 = np.mean([TGN_Jun['K1amp'],TGS_Feb['K1amp'],TGS_Jun['K1amp']]) * K1_amp_finecorrection
elev_pha_K1 = np.mean([TGN_Jun['K1pha_GMT'],TGS_Feb['K1pha_GMT'],
TGS_Jun['K1pha_GMT']]) + K1_pha_finecorrection
Z1_K1 = elev_amp_K1*np.cos(np.radians(elev_pha_K1))
Z2_K1 = elev_amp_K1*np.sin(np.radians(elev_pha_K1))
#M2 surface elevation
elev_amp_M2 = np.mean([TGN_Jun['M2amp'],TGS_Feb['M2amp'],TGS_Jun['M2amp']])*M2_amp_finecorrection
elev_pha_M2 = np.mean([TGN_Jun['M2pha_GMT'],TGS_Feb['M2pha_GMT'],
TGS_Jun['M2pha_GMT']]) + M2_pha_finecorrection
print (elev_pha_K1)
Z1_M2 = elev_amp_M2*np.cos(np.radians(elev_pha_M2))
Z2_M2 = elev_amp_M2*np.sin(np.radians(elev_pha_M2))
Start to make the netcdf file
In [26]:
#reload(webtide_tools)
#wait.... what about the U and V being on different grids?
#.... nothing really changes because we don't have great spatial resolution of measurements...
#... perhaps we'll need to make an extra cell for the V grid....
#currents
webtide_tools.create_northern_tides(UZ1_K1, UZ2_K1,'U','K1', tidesetname, name='SalishSea')
webtide_tools.create_northern_tides(VZ1_K1, VZ2_K1,'V','K1', tidesetname, name='SalishSea')
webtide_tools.create_northern_tides(UZ1_M2, UZ2_M2,'U','M2', tidesetname, name='SalishSea')
webtide_tools.create_northern_tides(VZ1_M2, VZ2_M2,'V','M2', tidesetname, name='SalishSea')
#surface elevation
webtide_tools.create_northern_tides(Z1_K1, Z2_K1,'T','K1', tidesetname, name='SalishSea')
webtide_tools.create_northern_tides(Z1_M2, Z2_M2,'T','M2', tidesetname, name='SalishSea')
Now, remake the tides in the johnstone_strait_contd notebook.
In [ ]: