Notebook

Comparing merge-tests on Salish and Japser. A week long run to see how the differences evolve.

In [1]:

%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import matplotlib as mpl
import netCDF4 as NC
import numpy as np
import scipy.interpolate as sp
import math

In [2]:

def combine_data(data_list):
    #combines ouput files so that they are easier to analyze
    #data_list is a dict object that contains the netcdf handles for the files of importance. 
    us={}; vs={}; lats={}; lons={}; sals={}; tmps={}; sshs={};
    for k in data_list:
        net=data_list.get(k)
        us[k]=net.variables['vozocrtx']
        vs[k]=net.variables['vomecrty']
        lats[k]=net.variables['nav_lat']        
        lons[k]=net.variables['nav_lon']
        tmps[k]=net.variables['votemper']
        sals[k]=net.variables['vosaline']
        sshs[k]=net.variables['sossheig']
    return us, vs, lats, lons, tmps, sals, sshs

In [3]:

#load up the bathymetry. 
grid = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
bathy = grid.variables['Bathymetry'][:,:]
X = grid.variables['nav_lon'][:,:]
Y = grid.variables['nav_lat'][:,:]

In [4]:

#load in the datas. Trying something new with dict objects...
r1 = 'r3942-4510_long';
runname1 = '/data/nsoontie/MEOPAR/SalishSea/results/merge-tests/' + r1 
print runname1

run1_thal={}
for k in range(1,7):
    string = 'thal' + str(k)
    run1_thal[string] = NC.Dataset(runname1 +'/1h_Thalweg' +str(k)+ '.nc','r');

#second dataset
r2 = 'r3942-4510Jasper_long'
runname2 = '/data/nsoontie/MEOPAR/SalishSea/results/merge-tests/' +r2
print runname2

run2_thal={}
for k in range(1,7):
    string = 'thal' + str(k)
    run2_thal[string] = NC.Dataset(runname2 +'/1h_Thalweg' +str(k)+ '.nc','r');

[us1, vs1, lats1, lons1, tmps1, sals1, sshs1] = combine_data(run1_thal)
[us2, vs2, lats2, lons2, tmps2, sals2, sshs2] = combine_data(run2_thal)

#finally grab a times placer and depths
net=run1_thal.get('thal1')
times = net.variables['time_counter']
depths = net.variables['deptht']

/data/nsoontie/MEOPAR/SalishSea/results/merge-tests/r3942-4510_long
/data/nsoontie/MEOPAR/SalishSea/results/merge-tests/r3942-4510Jasper_long

Before plotting anything, lets look at the points on the bathymetry.

In [5]:

#Plot bathy
ax=[-126.3,-122.1,47,51]

istorms=np.array([329, 196, 214, 124]);
jstorms=np.array([469, 299, 352, 748]);
iplume = 290
jplume = 440

fig =plt.figure(figsize=(6,8))
pylab.pcolor(X,Y,bathy)
pylab.colorbar()
for key in lats1:
    pylab.plot(lons1[key][0,0], lats1[key][0,0],'o',label=key)
pylab.plot(X[jstorms-1,istorms-1],Y[jstorms-1,istorms-1], 'ko')
pylab.plot(X[jplume-1,iplume-1],Y[jplume-1,iplume-1], 'ko')
pylab.axis(ax)
pylab.legend(loc=0)

Out[5]:

<matplotlib.legend.Legend at 0x50186d0>

Now lets try to look at some of the data. Lets plot the difference in sea surface height over time at each thalweg station

In [6]:

plt.figure(figsize=(20,8))
panel=1
ax=[0,7,-0.002,0.002];
for key in run1_thal:
    pylab.subplot(3,2,panel)
    ssh_diff = sshs1[key][:,0,0] - sshs2[key][:,0,0];
    #pylab.plot(times[:]/86400,sshs1[key][:,0,0],label=r1)
    #pylab.plot(times[:]/86400,sshs2[key][:,0,0],label=r2)
    pylab.plot(times[:]/86400,ssh_diff,label=r1 +' - ' +r2)
    pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')')
    pylab.ylabel('ssh (m)')
    pylab.axis(ax)
    panel = panel+1
pylab.legend(loc=2)    

Out[6]:

<matplotlib.legend.Legend at 0x2885b10>

There are pretty small (less than 1 cm).

Same for salinities now. Try to fit in surface, 25m depth and 100 m depth all in the same plot.

In [7]:

plt.figure(figsize=(20,8))
panel=1
ax=[0,7,-0.05,0.05];
depthlevel=0;
print 'Plotting diffefrence in salinity: ' +r1 +' - ' +r2
for key in run1_thal:
    pylab.subplot(3,2,panel)
    sal_diff = sals1[key][:,depthlevel,0,0] - sals2[key][:,depthlevel,0,0];
    #pylab.plot(times[:]/86400,sals1[key][:,depthlevel,0,0],label=r1)
    #pylab.plot(times[:]/86400,sals2[key][:,depthlevel,0,0],label=r2)
    pylab.plot(times[:]/86400,sal_diff,label = 'surface')
    pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')')
    pylab.ylabel('Salinity')
    pylab.axis(ax)
    panel = panel+1
    
# do it again for 25m depth
panel = 1
depthlevel=25
for key in run1_thal:
    pylab.subplot(3,2,panel)
    sals1P=np.zeros(times.shape[0]); sals2P=np.zeros(times.shape[0]); 
    # interpolate to 25m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],sals1[key][t,:,0,0]); sals1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],sals2[key][t,:,0,0]); sals2P[t] = f(depthlevel)
    sal_diff = sals1P-sals2P
    #pylab.plot(times[:]/86400,sals1P,'b*',label=depthlevel)
    #pylab.plot(times[:]/86400,sals2P,'g*')
    pylab.plot(times[:]/86400,sal_diff, label = depthlevel)
    panel = panel+1
    
# do it again for 100m depth
panel = 1
depthlevel=100
for key in run1_thal:
    pylab.subplot(3,2,panel)
    sals1P=np.zeros(times.shape[0]); sals2P=np.zeros(times.shape[0]); 
    # interpolate to 25m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],sals1[key][t,:,0,0]); sals1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],sals2[key][t,:,0,0]); sals2P[t] = f(depthlevel)
    sal_diff = sals1P-sals2P
    #pylab.plot(times[:]/86400,sals1P,'bs',label=depthlevel)
    #pylab.plot(times[:]/86400,sals2P,'gs')
    pylab.plot(times[:]/86400,sal_diff,label = depthlevel)
    panel = panel+1
    
pylab.legend(loc=2)

Plotting diffefrence in salinity: r3942-4510_long - r3942-4510Jasper_long

Out[7]:

<matplotlib.legend.Legend at 0xdf26d90>

Again, this is pretty small compared to the actual salinity values. There are larger differences in thalweg 3 (mixing area). Also, in thalweg 1 (close to JdF).

Let's look at u now.

In [8]:

plt.figure(figsize=(20,8))
panel=1
ax=[0,7,-.005,.005];
depthlevel=0;
for key in run1_thal:
    pylab.subplot(3,2,panel)
    u_diff = us1[key][:,depthlevel,0,0] - us2[key][:,depthlevel,0,0];
    #pylab.plot(times[:]/86400,us1[key][:,depthlevel,0,0],label=r1)
    #pylab.plot(times[:]/86400,us2[key][:,depthlevel,0,0],label=r2)
    pylab.plot(times[:]/86400,u_diff,label='surface')
    pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')')
    pylab.ylabel('u (m/s)')
    pylab.axis(ax)
    panel = panel+1
    
# do it again for 25m depth
panel = 1
depthlevel=25
for key in run1_thal:
    pylab.subplot(3,2,panel)
    us1P=np.zeros(times.shape[0]); us2P=np.zeros(times.shape[0]); 
    # interpolate to 25m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],us1[key][t,:,0,0]); us1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],us2[key][t,:,0,0]); us2P[t] = f(depthlevel)
    u_diff = us1P-us2P
    #pylab.plot(times[:]/86400,us1P,'b*',label=depthlevel)
    #pylab.plot(times[:]/86400,us2P,'g*')
    pylab.plot(times[:]/86400,u_diff,label=depthlevel)
    panel = panel+1
    
# do it again for 100m depth
panel = 1
depthlevel=100
for key in run1_thal:
    pylab.subplot(3,2,panel)
    us1P=np.zeros(times.shape[0]); us2P=np.zeros(times.shape[0]); 
    # interpolate to 25m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],us1[key][t,:,0,0]); us1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],us2[key][t,:,0,0]); us2P[t] = f(depthlevel)
    u_diff = us1P-us2P
    #pylab.plot(times[:]/86400,us1P,'bs',label=depthlevel)
    #pylab.plot(times[:]/86400,us2P,'gs')
    pylab.plot(times[:]/86400,u_diff,label=depthlevel)
    panel = panel+1
    
pylab.legend(loc=2)    

Out[8]:

<matplotlib.legend.Legend at 0xed2c790>

At these thalweg stations, once again the differences are pretty small (less than 1 cm/s). The largest differences are in thalweg 1 and thalweg 3.

Ok, now for v.

In [9]:

plt.figure(figsize=(20,8))
panel=1
ax=[0,7,-0.005,0.005];
depthlevel=0;
for key in run1_thal:
    pylab.subplot(3,2,panel)
    v_diff = vs1[key][:,depthlevel,0,0] - vs2[key][:,depthlevel,0,0];
    #pylab.plot(times[:]/86400,vs1[key][:,depthlevel,0,0],label=r1)
    #pylab.plot(times[:]/86400,vs2[key][:,depthlevel,0,0],label=r2)
    pylab.plot(times[:]/86400,v_diff,label='surface')
    pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')')
    pylab.ylabel('v (m/s)')
    pylab.axis(ax)
    panel = panel+1
    
# do it again for 25m depth
panel = 1
depthlevel=25
for key in run1_thal:
    pylab.subplot(3,2,panel)
    vs1P=np.zeros(times.shape[0]); vs2P=np.zeros(times.shape[0]); 
    # interpolate to 25m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],vs1[key][t,:,0,0]); vs1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],vs2[key][t,:,0,0]); vs2P[t] = f(depthlevel)
    v_diff = vs1P-vs2P
    #pylab.plot(times[:]/86400,vs1P,'b*',label=depthlevel)
    #pylab.plot(times[:]/86400,vs2P,'g*')
    pylab.plot(times[:]/86400,v_diff,label=depthlevel)
    panel = panel+1
    
# do it again for 100m depth
panel = 1
depthlevel=100
for key in run1_thal:
    pylab.subplot(3,2,panel)
    vs1P=np.zeros(times.shape[0]); vs2P=np.zeros(times.shape[0]); 
    # interpolate to 100m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],vs1[key][t,:,0,0]); vs1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],vs2[key][t,:,0,0]); vs2P[t] = f(depthlevel)
    v_diff = vs1P-vs2P
    #pylab.plot(times[:]/86400,vs1P,'bs',label=depthlevel)
    #pylab.plot(times[:]/86400,vs2P,'gs')
    pylab.plot(times[:]/86400,v_diff,label=depthlevel)
    panel = panel+1
    
pylab.legend(loc=2)

Out[9]:

<matplotlib.legend.Legend at 0xfa3c390>

One last time with temperature.

In [10]:

plt.figure(figsize=(20,8))
panel=1
ax=[0,7,-.01,0.01];
depthlevel=0;
for key in run1_thal:
    pylab.subplot(3,2,panel)
    T_diff = tmps1[key][:,depthlevel,0,0] - tmps2[key][:,depthlevel,0,0];
    pylab.plot(times[:]/86400,T_diff,label='surface')
    pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')')
    pylab.ylabel('Temperature')
    pylab.axis(ax)
    panel = panel+1
    
# do it again for 25m depth
panel = 1
depthlevel=25
for key in run1_thal:
    pylab.subplot(3,2,panel)
    tmps1P=np.zeros(times.shape[0]); tmps2P=np.zeros(times.shape[0]); 
    # interpolate to 25m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],tmps1[key][t,:,0,0]); tmps1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],tmps2[key][t,:,0,0]); tmps2P[t] = f(depthlevel)
    T_diff = tmps1P-tmps2P
    pylab.plot(times[:]/86400,T_diff,label=depthlevel)
    panel = panel+1
    
# do it again for 100m depth
panel = 1
depthlevel=100
for key in run1_thal:
    pylab.subplot(3,2,panel)
    tmps1P=np.zeros(times.shape[0]); tmps2P=np.zeros(times.shape[0]); 
    # interpolate to 25m first
    for t in range(0,times.shape[0]):
        f = sp.interp1d(depths[:],tmps1[key][t,:,0,0]); tmps1P[t] = f(depthlevel)
        f = sp.interp1d(depths[:],tmps2[key][t,:,0,0]); tmps2P[t] = f(depthlevel)
    T_diff = tmps1P-tmps2P
    pylab.plot(times[:]/86400,T_diff,label=depthlevel)
    panel = panel+1

pylab.legend(loc=2)

Out[10]:

<matplotlib.legend.Legend at 0x10746f50>

Look at the other points now. Let's just compare surface values.

In [11]:

stations = {'PointAtkinson': 1, 'Victoria': 1, 'PatriciaBay': 1, 'CampbellRiver': 1, 'Plume': 1}
run1_stations={}
run2_stations ={}


for key in stations:
    string = runname1 + '/1h_' + key + '.nc'
    run1_stations[key] = NC.Dataset(string,'r');
    string = runname2 + '/1h_' + key + '.nc'
    run2_stations[key] = NC.Dataset(string,'r');
    
[us1, vs1, lats1, lons1, tmps1, sals1, sshs1] = combine_data(run1_stations)
[us2, vs2, lats2, lons2, tmps2, sals2, sshs2] = combine_data(run2_stations)
print 'plotting ' + r1 +' - ' + r2

fig, ((ax_ssh, ax_sal), (ax_u, ax_v)) = plt.subplots(2, 2, figsize=(18, 10))
for locn in stations:
    ssh_diff=sshs1[locn][:,0,0] -sshs2[locn][:,0,0]
    u_diff=us1[locn][:,0,0,0] -us2[locn][:,0,0,0]
    v_diff=vs1[locn][:,0,0,0] -vs2[locn][:,0,0,0]
    sals_diff=sals1[locn][:,0,0,0] -sals2[locn][:,0,0,0]
    ax_ssh.plot(times[:]/86400, ssh_diff, marker='+', label=locn)
    ax_sal.plot(times[:]/86400,sals_diff, marker='+', label=locn)
    ax_u.plot(times[:]/86400,u_diff, marker='+', label=locn)
    ax_v.plot(times[:]/86400,v_diff, marker='+', label=locn)
for ax in (ax_ssh, ax_sal, ax_u, ax_v):
    start, end = ax.set_xlim(0, 7)
    ax.set_xlabel('time')
    

    
vmin=-0.005; vmax=0.005
ax_ssh.set_ylim(vmin,vmax); ax_ssh.set_ylabel('SSH (m)');
ax_u.set_ylim(-0.001,0.001); ax_u.set_ylabel('Surface u (m/s)');
ax_v.set_ylim(-0.001,0.001); ax_v.set_ylabel('Surface v (m/s)');
ax_sal.set_ylim(-.05,.05); ax_sal.set_ylabel('Surface Salinity');
ax_sal.legend(loc='best')

plotting r3942-4510_long - r3942-4510Jasper_long

Out[11]:

<matplotlib.legend.Legend at 0x119caf10>

Again, the differences are quite small. Let's look at the maximum differences over the whole domain. Here, I am looking at every grid cell and not just at the points along our stations.

In [12]:

fT1=  NC.Dataset(runname1 + '/SalishSea_4h_20020915_20020921_grid_T.nc','r');
fU1=  NC.Dataset(runname1 + '/SalishSea_4h_20020915_20020921_grid_U.nc','r');
fV1=  NC.Dataset(runname1 + '/SalishSea_4h_20020915_20020921_grid_V.nc','r');

fT2=  NC.Dataset(runname2 + '/SalishSea_4h_20020915_20020921_grid_T.nc','r');
fU2=  NC.Dataset(runname2 + '/SalishSea_4h_20020915_20020921_grid_U.nc','r');
fV2=  NC.Dataset(runname2 + '/SalishSea_4h_20020915_20020921_grid_V.nc','r');
    
u1 = fU1.variables['vozocrtx']; v1 = fV1.variables['vomecrty']; ssh1 = fT1.variables['sossheig'];

u2 = fU2.variables['vozocrtx']; v2 = fV2.variables['vomecrty']; ssh2 = fT2.variables['sossheig'];

times_full = fT2.variables['time_counter'];

In [13]:

#Calculate max difference over the field for all times.

fig =plt.figure(figsize=(12,12))
print ' Plotting difference in ssh over time: ' + r1 + ' - ' + r2

sshmaxdiff=np.zeros((42,1));    
for ii in range(0,42):
    sshdiff =ssh1[ii,:,:]-ssh2[ii,:,:];
    sshmaxdiff[ii] = abs(sshdiff[:]).max()
  
fig =plt.figure(figsize=(12,6))
pylab.plot(times_full[:]/86400,sshmaxdiff)
pylab.xlabel('time (d)')
pylab.ylabel('difference in ssh (m)')
pylab.title('Maximum difference in ssh')
    

 Plotting difference in ssh over time: r3942-4510_long - r3942-4510Jasper_long

Out[13]:

<matplotlib.text.Text at 0x1259fc90>

<matplotlib.figure.Figure at 0xd1c4690>

When I take into account the entire domain, the differences are more substantial. Here, the max difference in sea surface height is up to 2 cm. Over a longer time the behaviour seems to mellow out.

It isn't clear here where the max differences are occuring. Previous runs show large differences near JdF and in the mixing region.

Let's look at the other fields too just to check.

In [14]:

#Calculate max difference over the field for all times.

fig =plt.figure(figsize=(12,12))
print ' Plotting difference in u over time: ' + r1 + ' - ' + r2

umaxdiff=np.zeros((42,1));    
for ii in range(0,42):
    udiff =u1[ii,:,:,:]-u2[ii,:,:,:];
    umaxdiff[ii] = abs(udiff[:]).max()
  
fig =plt.figure(figsize=(12,6))
pylab.plot(times_full[:]/86400,umaxdiff)
pylab.xlabel('time (d)')
pylab.ylabel('difference in u (m/s)')
pylab.title('Maximum difference in u')
    

 Plotting difference in u over time: r3942-4510_long - r3942-4510Jasper_long

Out[14]:

<matplotlib.text.Text at 0x124308d0>

<matplotlib.figure.Figure at 0x10764b10>

This worried me a little: seeing current differences up to 0.8m/s! But again, these differences seem to settle out over time. There is also a weird daily signal in the first three days. I'm not sure what could be causing that.

Do this again for v.

In [15]:

#Calculate max difference over the field for all times.

vmaxdiff=np.zeros((42,1));    
for ii in range(0,42):
    vdiff =v1[ii,:,:,:]-v2[ii,:,:,:];
    vmaxdiff[ii] = abs(vdiff[:]).max()
  
    
print ' Plotting difference in v over time: ' + r1 + ' - ' + r2
fig =plt.figure(figsize=(12,6))
pylab.plot(times_full[:]/86400,vmaxdiff)
pylab.xlabel('time (d)')
pylab.ylabel('difference in v (m/s)')
pylab.title('Maximum difference in v')
    

 Plotting difference in v over time: r3942-4510_long - r3942-4510Jasper_long

Out[15]:

<matplotlib.text.Text at 0x139e4350>

Summary¶

Differences at the stations (thalwegs and storm surge points) are not very large.
When I take into account the entire domain, the maximum differences are much larger. However, these differences generally decrease over time and by the end of the week are pretty small. There is also a weird daily signal for the first three days, most apparent in u.
The adjustment to forcing is slightly different on the two machines. This is probably to be expected. It was reassuring to see these differences decrease (rather than increase) over time.

In [ ]: