<jupyter_start><jupyter_text>Get the forcing from WebTide (http://www.bio.gc.ca/science/research-recherche/ocean/webtide/index-eng.php) for the Juan de Fuca boundary of the model and created the web forcing files for the model, including all tidal corrections<jupyter_code>#preamble
%matplotlib inline
from matplotlib import pylab
import netCDF4 as NC
import csv 
import numpy
import importlib
from itertools import islice
from math import cos, sin, tan, radians
import sys
import webtide_tools
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}}}
# These are the corrections to be used.
tidesetname = 'N36_AF'
#list of constituents in Webtide predictions
constlist = ['O1', 'Q1', 'K1', 'P1', 'N2', 'M2', 'S2', 'K2']
# these are the Fine Tuning Correction Factors.
# A means All, A2 : all semi-diurnals, A1: all diurnals
CFactors = {'A2 Flux':1.2240, 'A2 Phase': TSin['West']['M2']['Pha'], 
            'A2 Amp': TSin['West']['M2']['Amp'], 'A2 Shift' : 0.,
            'S2 Phase': TSin['West']['S2']['Pha'], 
            'S2 Amp': TSin['West']['S2']['Amp'],
            'N2 Phase': TSin['West']['N2']['Pha'], 
            'N2 Amp': TSin['West']['N2']['Amp'],
            'K2 Phase': 0., 'K2 Amp': 0.,
            'A1 Phase': TSin['West']['K1']['Pha'], 
            'A1 Amp': TSin['West']['K1']['Amp'],
            'O1 Phase': TSin['West']['O1']['Pha'], 
            'O1 Amp': TSin['West']['O1']['Amp'],
            'Q1 Phase': TSin['West']['Q1']['Pha'], 
            'Q1 Amp': TSin['West']['Q1']['Amp']}
CFactors['K2 Phase'] = -5.7
CFactors['K2 Amp'] = 0.235
CFactors['P1 Phase'] = - 3.5
CFactors['P1 Amp'] = 0.31<jupyter_output><empty_output><jupyter_text>Get the co-ordinates of the boundary and put them in a text file for WebTide. Note that surface elevation, u and v are offset from each other on a Arakawa C grid. So the lat/lon of each grid will be different for each of them. We therefore need a different list of co-ordinates to give WebTide for T, U and V.Use new bathymetry with smoothed JdF mouth.K1 phase decreased by 5 degrees in U,V,eta.K1 amplitude increased by factor of 1.15 in U,V,etaM2 phase decreased by 5 degrees in U,V,eta<jupyter_code>#load in the bathymetry netcdf data: use old bathymetry as this is redoing unchanged calculations
fB = NC.Dataset('../../../NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')

#get the depth data from the netcdf
depth = fB.variables['Bathymetry'][:]
# masking just causes problems, set depth to 0 on land
depth[depth.mask] = 0
a, b = depth.shape
print(a,b)
#load in the co-ordinates file
fC = NC.Dataset('../../../NEMO-forcing/grid/coordinates_seagrid_SalishSea.nc','r')

#glam = longitude, gphi = latitude
lonT = fC.variables['glamt'][0,:,0]
latT = fC.variables['gphit'][0,:,0]
lonU = fC.variables['glamu'][0,:,0]
latU = fC.variables['gphiu'][0,:,0]
lonV = fC.variables['glamv'][0,:,0]
latV = fC.variables['gphiv'][0,:,0]

#check it's the right boundary :)
pylab.contourf(fC.variables['glamt'][0,:,:],fC.variables['gphit'][0,:,:],depth)
pylab.plot(lonT,latT,'oy')
#i and j indices of the western boundary
i = a-1
j = 0

#find the boundary from the depth data
depthT = depth[:,j]
depthU = 0.5*(depth[:,j]+depth[:,j+1])
depthV = numpy.zeros(len(depthT))
depthV[0:i-1] = 0.5*(depth[0:i-1,j]+depth[1:i,j])
depthV[i] = depth[i,j]

#plot all the depth data
#if the three depth variables line up then I've done it right!
pylab.plot(depthT,'x',depthU,'o',depthV,'*')
pylab.xlabel('grid number')
pylab.ylabel('depth [m]')
pylab.legend(('T','U','V'))
pylab.title('Check the three grids line up')

#find where lats and lons the depths are not 0 - this is the boundary
boundlatT = latT[depthT!=0]
boundlonT = lonT[depthT!=0]
boundlatU = latU[depthU!=0]
boundlonU = lonU[depthU!=0]
boundlatV = latV[depthV!=0]
boundlonV = lonV[depthV!=0]

print('number of T cells along boundary = '+str(len(boundlatT)))
print('number of U cells along boundary = '+str(len(boundlatU)))
print('number of V cells along boundary = '+str(len(boundlatV)))
# skippable
print('T')
I = numpy.where(depthT!=0)
print(I[0][-1]+6)
print(I[0][0]-5)

print(len(depth[I]))
print('u')
print(numpy.where(depthU!=0))
print('v')
print(numpy.where(depthV!=0))<jupyter_output>T
471
384
77
u
(array([389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401,
       402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414,
       415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427,
       428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440,
       441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453,
       454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465]),)
v
(array([388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400,
       401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413,
       414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426,
       427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439,
       440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452,
       453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465]),)<jupyter_text>skippableWrite the lats and lons of each grid to a text fileT gridcreate a .csv file of lat/lons for tidal elevation predictions, to import into Webtidewith open('jdf_boundary_coords_T.txt', 'wb') as csvfile:    writer = csv.writer(csvfile, delimiter=' ')    for t in numpy.arange(0,len(boundlatT)):            writer.writerow([boundlonT[t],boundlatT[t]])                     U gridcreate a .csv file of lat/lons for tidal elevation predictions, to import into Webtidewith open('jdf_boundary_coords_U.txt', 'wb') as csvfile:    writer = csv.writer(csvfile, delimiter=' ')    for t in numpy.arange(0,len(boundlatU)):            writer.writerow([boundlonU[t],boundlatU[t]])         V gridcreate a .csv file of lat/lons for tidal elevation predictions, to import into Webtidewith open('jdf_boundary_coords_V.txt', 'wb') as csvfile:    writer = csv.writer(csvfile, delimiter=' ')    for t in numpy.arange(0,len(boundlatV)):            writer.writerow([boundlonV[t],boundlatV[t]])          Now I start the WebTide application, select 'Import Tide Markers' then select 'Get Constituents'. An html page opens with the constituents for each of the points. I select all the data and copy into a .csv file. I do this separately for T, U and V grids. For the T grid, I get constituents for elevation. For the U grid, I get constituents for currents but I only use the u component. For the V grid, I get constituents for currents but I only use the v component.NOTE: I did this once for the subdomain grid, but the lat/lon of the boundary points didn't change so I didn't change it for the bigger domain.  Read in Webtide data, change from amplitude/phase to Z1/Z2 (required by NEMO) and create NetCDF files<jupyter_code># skippable

importlib.reload(webtide_tools)
# Use the flux/amplitude/phase corrections and calculate and write the files.


for i in numpy.arange(0,len(constlist)):
    
    # get the predictsion results for elevation, y and v and convert to required format
    # the last number, puts a number after SalishSea.  Here a "2" says the bathymetry with
    # smoothed JdF mouth
    # if its longer than a single digit the whole name is used instead of SalishSea
    elZ1, elZ2 = webtide_tools.create_tide_netcdf('T', constlist[i], depthT, 'SalishSea', tidesetname, CFactors)
    uZ1, uZ2 = webtide_tools.create_tide_netcdf('U', constlist[i], depthU, 'SalishSea', tidesetname, CFactors)
    vZ1, vZ2 = webtide_tools.create_tide_netcdf('V', constlist[i], depthV, 'SalishSea', tidesetname, CFactors)
    
    #plot the results for each constituent to check it makes sense
    pylab.figure()
    #define some arbitrary time to check results
    theta1 = numpy.pi/4

    #plot the elevation results
    pylab.subplot(4,1,1)
    Ael = numpy.sqrt(elZ1**2+elZ2**2)
    phiel = numpy.arctan2(elZ2,elZ1)
    elN = Ael*numpy.cos(theta1-phiel)
    pylab.title(constlist[i])
    pylab.plot(elN,'r')
    pylab.ylabel('elevation (m)')
    pylab.ylim((-0.4,0.4))
    
    #calculate the actual components u and v
    Au = numpy.sqrt(uZ1**2+uZ2**2)
    Av = numpy.sqrt(vZ1**2+vZ2**2)
    phiu = numpy.arctan2(uZ2,uZ1)
    phiv = numpy.arctan2(vZ2,vZ1)
    uN = Au*numpy.cos(theta1-phiu)
    vN = Av*numpy.cos(theta1-phiv)
    
    #plot the u and v velocities separately
    pylab.subplot(4,1,2)
    pylab.plot(uN,'b')
    pylab.ylabel('u (m/s)')
    pylab.ylim((-0.3,0.3))

    pylab.subplot(4,1,3)
    pylab.plot(vN,'m')
    pylab.ylabel('v (m/s)')
    pylab.ylim((-0.3,0.3))
    
    #make a little quiver plot to check this makes sense
    pylab.subplot(4,1,4)
    pylab.quiver([0]*78,numpy.linspace(0,77,78),uN[0:76],vN[0:76],scale=1)
    
    if constlist[i] == "M2":
        phase_shift = (phiu-phiel)*180./numpy.pi
        phase_u = (phiu)*180./numpy.pi
        phase_T = (phiel)*180./numpy.pi

#    print phiel*180/numpy.pi
# skippable
pylab.plot(phase_u)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1591542031975245, 'fraction_numerical': 0.11397627643115008, 'mean_word_length': 4.095112979506043, 'pattern_counts': {'":': 1, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '198780', 'license': 'ECL-2.0', 'n_tokens_mistral': 4532, 'n_tokens_neox': 3714, 'n_words': 1063, 'path': 'I_ForcingFiles/Tides/webtide_forcing.ipynb', 'repo_name': 'mdunphy/tools', 'size': '9695'}