A code to look at the data from the crashed run. What is going wrong! Looks at currents, bathymetry, and sea surface height around the crashed time step.
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.
f = NC.Dataset('/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/nov2009/all_forcing_crash/output.abort.nc','r');
#fU2 = NC.Dataset('/data/nsoontie/MEOPAR/SalishSea/results/viscosity/nu40_crash/SalishSea_4h_20020915_20020921_grid_U.nc','r');
grid = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea.nc','r')
bathy = grid.variables['Bathymetry'][:,:]
In [3]:
#Print some info about the variables.
for dim in f.dimensions.values():
print dim
<type 'netCDF4.Dimension'>: name = 'x', size = 398 <type 'netCDF4.Dimension'>: name = 'y', size = 898 <type 'netCDF4.Dimension'>: name = 'deptht', size = 40 <type 'netCDF4.Dimension'> (unlimited): name = 'time_counter', size = 1
In [4]:
f.variables.keys()
Out[4]:
[u'nav_lon', u'nav_lat', u'deptht', u'time_counter', u'vosaline', u'votemper', u'sossheig', u'vozocrtx', u'vomecrty', u'vovecrtz', u'sowaflup', u'sohefldo', u'soshfldo', u'soicecov', u'sozotaux', u'sometauy']
In [5]:
for var in f.variables.values():
print var
<type 'netCDF4.Variable'>
float32 nav_lon(y, x)
standard_name: longitude
units: degrees_east
valid_min: -126.4
valid_max: -121.318
long_name: Longitude
nav_model: Default grid
unlimited dimensions:
current shape = (898, 398)
filling off
<type 'netCDF4.Variable'>
float32 nav_lat(y, x)
standard_name: latitude
units: degrees_north
valid_min: 46.8597
valid_max: 51.1048
long_name: Latitude
nav_model: Default grid
unlimited dimensions:
current shape = (898, 398)
filling off
<type 'netCDF4.Variable'>
float32 deptht(deptht)
axis: Z
standard_name: model_level_number
units: m
positive: down
valid_min: 0.5
valid_max: 441.466
title: deptht
long_name: Vertical T levels
unlimited dimensions:
current shape = (40,)
filling off
<type 'netCDF4.Variable'>
float64 time_counter(time_counter)
axis: T
standard_name: time
units: seconds since 2009-11-14 22:13:10
calendar: gregorian
title: Time
long_name: Time axis
time_origin: 2009-NOV-14 22:13:10
unlimited dimensions: time_counter
current shape = (1,)
filling off
<type 'netCDF4.Variable'>
float32 vosaline(time_counter, deptht, y, x)
units: PSU
standard_name: Salinity
_FillValue: 9.96921e+36
long_name: Salinity
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter deptht nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 40, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 votemper(time_counter, deptht, y, x)
units: C
standard_name: Temperature
_FillValue: 9.96921e+36
long_name: Temperature
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter deptht nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 40, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 sossheig(time_counter, y, x)
units: m
standard_name: Sea Surface Height
_FillValue: 9.96921e+36
long_name: Sea Surface Height
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 vozocrtx(time_counter, deptht, y, x)
units: m/s
standard_name: Zonal Current
_FillValue: 9.96921e+36
long_name: Zonal Current
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter deptht nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 40, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 vomecrty(time_counter, deptht, y, x)
units: m/s
standard_name: Meridional Current
_FillValue: 9.96921e+36
long_name: Meridional Current
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter deptht nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 40, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 vovecrtz(time_counter, deptht, y, x)
units: m/s
standard_name: Vertical Velocity
_FillValue: 9.96921e+36
long_name: Vertical Velocity
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter deptht nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 40, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 sowaflup(time_counter, y, x)
units: Kg/m2/S
standard_name: Net Upward Water Flux
_FillValue: 9.96921e+36
long_name: Net Upward Water Flux
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 sohefldo(time_counter, y, x)
units: W/m2
standard_name: Net Downward Heat Flux
_FillValue: 9.96921e+36
long_name: Net Downward Heat Flux
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 soshfldo(time_counter, y, x)
units: W/m2
standard_name: Shortwave Radiation
_FillValue: 9.96921e+36
long_name: Shortwave Radiation
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 soicecov(time_counter, y, x)
units: [0,1]
standard_name: Ice fraction
_FillValue: 9.96921e+36
long_name: Ice fraction
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 sozotaux(time_counter, y, x)
units: N/m2
standard_name: Zonal Wind Stress
_FillValue: 9.96921e+36
long_name: Zonal Wind Stress
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 898, 398)
filling off
<type 'netCDF4.Variable'>
float32 sometauy(time_counter, y, x)
units: N/m2
standard_name: Meridional Wind Stress
_FillValue: 9.96921e+36
long_name: Meridional Wind Stress
online_operation: inst(x)
interval_operation: 10.0
interval_write: 10.0
coordinates: time_counter nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (1, 898, 398)
filling off
In [6]:
u_vel = f.variables['vozocrtx'] #u currents
t=f.variables['time_counter']
ssh=f.variables['sossheig']
sal=f.variables['vosaline']
v_vel=f.variables['vomecrty'] # v currents
w_vel=f.variables['vovecrtz']
#bathy lat lon
b_lat = grid.variables['nav_lat']
b_lon = grid.variables['nav_lon']
In [7]:
#quick plot of the u field
ibad=220; jbad=184
depth=3
plt.figure(figsize=(14,8))
cmin=-2; cmax=2;
ax=[200,250,150,220]
#print time
tnow=t[0]/86400
print tnow
#plot U
U=u_vel[0,depth,:,:]
print U.min()
print U.max()
mu =U == 0
U= np.ma.array(U,mask=mu)
pylab.subplot(1,2,1)
pylab.pcolor(U,vmin=cmin,vmax=cmax)
#pylab.plot(ibad-1,jbad-1,marker='o') #marker for unstable point in the ocean.output file
pylab.colorbar()
#pylab.axis(ax)
#plot V
pylab.subplot(1,2,2)
V=v_vel[0,depth,:,:]
mu =V == 0
V= np.ma.array(V,mask=mu)
pylab.pcolor(V,vmin=cmin,vmax=cmax)
pylab.colorbar()
pylab.plot(ibad-1,jbad-1,marker='o') #marker for unstable point in the ocean.output file
#pylab.axis(ax)
print u_vel[0,3,jbad-1,ibad-1]
print v_vel[0,3,jbad-1,ibad-1]
0.0743055555556 -5.2533 25.6648 25.6648 0.0
In [8]:
#Vertical velocity
cmin=-0.005; cmax=0.005
depth=3
W=w_vel[0,depth,:,:]
print W.min()
print W.max()
mu =W == 0
W= np.ma.array(W,mask=mu)
pylab.pcolor(W,vmin=cmin,vmax=cmax)
pylab.colorbar()
pylab.axis(ax)
print w_vel[0,3,jbad-1,ibad-1]
-0.00915483 0.0180596 6.05582e-09
In [9]:
# load the bathymetry to look for shallow depths in this region.
plt.figure(figsize=(14,8))
cmin=0; cmax=5;
mu =bathy == 0
bathy= np.ma.array(bathy,mask=mu)
pylab.pcolor(bathy,vmin=cmin,vmax=cmax)
pylab.colorbar()
pylab.plot(ibad-1,jbad-1,marker='o') #marker for unstable point in the ocean.output file
pylab.axis(ax)
Out[9]:
[200, 250, 150, 220]
In [10]:
# look at sea surface height
plt.figure(figsize=(14,8))
mu =ssh == 0
ssh= np.ma.array(ssh,mask=mu)
pylab.pcolor(ssh[0,:,:],vmin=-4,vmax=4)
pylab.colorbar()
pylab.plot(ibad-1,jbad-1,marker='o') #marker for unstable point in the ocean.output file
#pylab.axis(ax)
print ssh.min()
print ssh.max()
print bathy[jbad-1,ibad-1]
print ssh[0,jbad-1,ibad-1]
-4.66426 1.90734 4.0 -3.2886
ssh is pretty negative here
In [11]:
# look at salinity
plt.figure(figsize=(14,8))
mu =sal == 0
sal= np.ma.array(sal,mask=mu)
pylab.pcolor(sal[0,depth,:,:],vmin=0,vmax=32)
pylab.colorbar()
pylab.plot(ibad-1,jbad-1,marker='o') #marker for unstable point in the ocean.output file
#pylab.axis(ax)
print sal[0,3,jbad-1,ibad-1]
29.3682
In [ ]: