An analysis of the open boundary conditions at the north.

These runs have been crashing, so look at an abort file:

In [1]:

%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import netCDF4 as NC
import numpy as np

In [2]:

# Load the data. Path name can be changed to look at different data.  
path = '/data/nsoontie/MEOPAR/SalishSea/results/north/tides/'

fU = NC.Dataset(path +'SalishSea_4h_20020915_20020915_grid_U.nc','r');
fV = NC.Dataset(path +'SalishSea_4h_20020915_20020915_grid_V.nc','r');
fT = NC.Dataset(path +'SalishSea_4h_20020915_20020915_grid_T.nc','r');
grid = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
bathy = grid.variables['Bathymetry'][:,:]

In [3]:

u_vel = fU.variables['vozocrtx']  #u currents
t=fU.variables['time_counter']
ssh=fT.variables['sossheig']
v_vel=fV.variables['vomecrty'] # v currents
sal = fT.variables['vosaline']
depths = fT.variables['deptht']
temp=fT.variables['votemper']

First, let's look at the bathymetry near the Johnstone Strait boundary. NEMO suggests that the depth should be constant over the first 4 grid cells. We may consider some smoothing as we did with JdF.

In [19]:

plt.figure(figsize=(8,8))
ax=[30, 100,880,898]
cmin=0; cmax=400;
mu =bathy == 0
bathy= np.ma.array(bathy,mask=mu)
pylab.pcolor(bathy,vmin=cmin,vmax=cmax)
pylab.colorbar()
pylab.axis(ax)

Out[19]:

[30, 100, 880, 898]

Now let's look at the u and v velocity fields, concentrating at the North in the lower plots.

In [28]:

#quick plot of the u field
depth=0
ts = 4

plt.figure(figsize=(10,10))
cmin=-2; cmax=2;

#plot U
U=u_vel[ts,depth,:,:]
print 'min/max U'
print U.min()
print U.max()
mu =U == 0
U= np.ma.array(U,mask=mu)
pylab.subplot(2,2,1)
pylab.pcolor(U,vmin=cmin,vmax=cmax)
pylab.colorbar()
pylab.title('U')
#pylab.axis(ax)

#plot V
pylab.subplot(2,2,2)
V=v_vel[ts,depth,:,:]
mu =V == 0
V= np.ma.array(V,mask=mu)
pylab.pcolor(V,vmin=cmin,vmax=cmax)
pylab.colorbar()
pylab.title('V')
#pylab.axis(ax)

print 'min/max V'
print V.min()
print V.max()

#plot U
mu =U == 0
U= np.ma.array(U,mask=mu)
pylab.subplot(2,2,3)
pylab.pcolor(U,vmin=-1,vmax=1)
pylab.colorbar()
pylab.title('U')
pylab.axis(ax)

#plot V
pylab.subplot(2,2,4)
V= np.ma.array(V,mask=mu)
pylab.pcolor(V,vmin=-1,vmax=1)
pylab.colorbar()
pylab.title('V')
pylab.axis(ax)

min/max U
-1.83501
2.08435
min/max V
-2.11239
2.38542

Out[28]:

[30, 100, 880, 898]

I want to check the tide forcing files to see if we have to make the correction that Susan made to her T+S. I don't think we accounted for Fortran counting when making the forcing files.

Quick look at the surface salinity and temperature.

In [29]:

smin=20; smax=32;
tmin=7; tmax=25;

plt.figure(figsize=(16,8))

S=sal[ts,depth,:,:]
mu =S == 0
S= np.ma.array(S,mask=mu)
pylab.subplot(1,2,1)
pylab.pcolor(S,vmin=smin,vmax=smax,cmap='winter')
pylab.colorbar()
#pylab.axis(ax)

pylab.subplot(1,2,2)
T=temp[ts,depth,:,:]
mu =T == 0
T= np.ma.array(T,mask=mu)
pylab.pcolor(T,vmin=tmin,vmax=tmax,cmap='hot')
pylab.colorbar()
#pylab.axis(ax)

Out[29]:

<matplotlib.colorbar.Colorbar instance at 0x4414d6c8>

SSH

In [40]:

emin=-2;emax=2;

plt.figure(figsize=(8,8))

E=ssh[ts,:,:]
mu =E == 0
E= np.ma.array(E,mask=mu)
pylab.pcolor(E,vmin=emin,vmax=emax,cmap='jet')
pylab.colorbar()
#pylab.axis(ax)

Out[40]:

<matplotlib.colorbar.Colorbar instance at 0x36868368>

I should look for a simulation without northern tides for comparison. I think the elevation in Johnstone is lower than other cases. Perhaps I should also reconsider the atmospheric forcing. I could turn off for now to isolate the effects of the tides.

SSH at North over time.

In [41]:

pylab.plot(ssh[:,895,45])
pylab.xlabel('output time')
pylab.ylabel('ssh (m)')
pylab.axis([0,5,-1.5,1.5])

Out[41]:

[0, 5, -1.5, 1.5]

Summary¶

This simulation used:

North - flather with tidal forcing (barotropic), relaxation to initial conditions (T+S), nothing (baroclinic)

West - flather with tides+external SSH clim (barotropic), relaxation to Masson clim (T+S), zero gradient (baroclinic)

viscosity = 100

Susan's suggestion for turning off the T+S boundary forcing and her corrections to the namelist worked nicely. This run was stable for one day. I think I am ready to lower viscosity to 50 and run for longer...

Note: it may be useful to output the hourly stations as well. At least to look at Campbell River.

Questions:

How can I tell if the forcing at the North is what I expect?
How does the forcing at the North affect the results? It may be worthwhile to design a control experiment.
Do I have to correct the tidal forcing files to account for Fortran counting like Susan did?

In [ ]: