Notebook to prepare NEMO initial stratification files from a single T,S profile from September 2003 taken at S4-1.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as NC
from scipy.interpolate import interp1d
from salishsea_tools import nc_tools
In [6]:
filedata = np.genfromtxt('../../../nemo-forcing/initial_strat/sg0318006.cnv', skiprows=98)
data = np.zeros((398+1, 12))
data[1:] = filedata
data[0] = data[1]
data[0,1] = 0.
pressure = data[:,1]

The twelve columns in data are: scan, pressure, temperature, conductivity, chl fluoresence, light transmission, PAR, oxygen, salinity, sigma, number of bins used to average and a flag.

Now we sort based on sigma (column 9, zero based) to get rid of any inversions and then put the pressure back in the correct order.

In [7]:
sorted = data[data[:,9].argsort()]
sorted[:,1] = pressure

Extract the other columns we are interested in

In [8]:
temp = sorted[:,2]
sal = sorted[:,8]
sigma = sorted[:,9]

Take a look at what we have done. Changes in Sigma (vertical axis) as a function of depth (horizontal axis)

In [9]:
plt.plot(sigma-data[:,9])
Out[9]:
[<matplotlib.lines.Line2D at 0x105eb03d0>]

Now we open the bathymetry file to give us latitude/longitude etc to base our input file upon. We could use the coordinates file, but its just more complicated. Use new topography with smoothed Juan de Fuca mouth.

In [10]:
fB = NC.Dataset('../../../nemo-forcing/grid/bathy_meter_SalishSea2.nc','r')
print fB.file_format
depth = fB.variables['Bathymetry']
lat = fB.variables['nav_lat']
lon = fB.variables['nav_lon']
BX = fB.dimensions['x']
print len(BX)
BY = fB.dimensions['y']
NETCDF4
398

Open a data file to determine the depths

In [22]:
fT = NC.Dataset('../../../nemo-forcing/grid/grid_bathy.nc')
nc_tools.show_dataset_attrs(fT)
deptht = fT.variables['deptht']
depth_out = deptht[:]
depths = fT.variables['grid_bathy']
depth = depths[:]
print depth.shape
file format: NETCDF4
Conventions: CF-1.6
title: NEMO z-partial-step Grid Level Depths
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://bitbucket.org/salishsea/tools/src/tip/bathymetry/NEMO-GridBathy.ipynb
references: https://bitbucket.org/salishsea/nemo-forcing/src/tip/grid/grid_bathy.nc
history: [2014-01-14 13:53:09] Created netCDF4 zlib=True dataset.
comment: Calculated by a 1-step NEMO run.
(40, 898, 398)

Interpolate Salinity and Temperature onto our Depths, note that below the maximum depth of the profile, it puts the maximum depth salinity and temperature. If you call it for depths above 0m, it will do the same. Don't do that!!!

In [16]:
print sal[-1]
salinity_function = interp1d(pressure, sal, bounds_error=False, fill_value=sal[-1])
temperature_function = interp1d(pressure, temp, bounds_error=False, fill_value=temp[-1])
31.3702

Set up our new NetCDF file (salinity first)

In [23]:
nemo = NC.Dataset('SS2_SoG0318_1y_salinity_nomask.nc', 'w', zlib=True)
# dataset attributes
nc_tools.init_dataset_attrs(
    nemo, 
    title='Salinity Initial Condition', 
    notebook_name='PrepareTS', 
    nc_filepath='../../../nemo-forcing/initial_strat/SS2_SoG0318_1y_salinity_nomask.nc',
    comment='Salinity from S4-1 September 2003, STRATOGEM Cruise 03-18; used at all grid points')
file format: NETCDF4
Conventions: CF-1.6
title: Salinity Initial Condition
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://bitbucket.org/salishsea/tools/src/tip/PrepareTS.ipynb
references: https://bitbucket.org/salishsea/nemo-forcing/src/tip/initial_strat/SS2_SoG0318_1y_salinity_nomask.nc
history: [2014-01-24 16:45:11] Created netCDF4 zlib=True dataset.
comment: Salinity from S4-1 September 2003, STRATOGEM Cruise 03-18; used at all grid points

Set up our dimensions. x and y match the bathymetry field, deptht and time counter is just 1

In [24]:
nemo.createDimension('x',len(BX))
nemo.createDimension('y',len(BY))
nemo.createDimension('deptht',size = len(depth_out))
nemo.createDimension('time_counter', None)
Out[24]:
<netCDF4.Dimension at 0x104b94c80>

Set up all our variables. Just copy latitude and longitude info.

In [25]:
# variables
nav_lat = nemo.createVariable('nav_lat', 'float32', ('y','x'))
nav_lat.long_name = 'Latitude'
nav_lat.units = 'degrees_north'
nav_lat = lat
nav_lon = nemo.createVariable('nav_lon', 'float32', ('y','x'))
nav_lon.long_name = 'Longitude'
nav_lon.units = 'degrees_east'
nav_lon = lon
deptht = nemo.createVariable('deptht', 'float32', ('deptht'))
deptht.long_name = 'Vertical T Levels'
deptht.units = 'm'
deptht.positive = 'down'
deptht.valid_range = np.array((4., 428.))
deptht = depth_out
time_counter = nemo.createVariable('time_counter', 'float32', ('time_counter'))
time_counter.units = 'seconds since 2003-09-10 12:27:00'
time_counter.long_name = 'Time axis'
time_counter.calendar = 'noleap'
vosaline = nemo.createVariable('vosaline', 'float32', 
                               ('time_counter','deptht','y','x'))
vosaline.units = 1
vosaline.long_name = 'Practical Salinity'  
vosaline.coordinates = 'nav_lon nav_lat deptht time_counter'
vosaline.grid = 'SalishSea2'

Assign time and salinity values:

In [28]:
print len(BX), len(BY)
vosaline[0,0,:,:] = temp[0]
for id in range(1,len(deptht)):
    vosaline[0,id,:,:] = salinity_function(deptht[id])
398 898
In [43]:
# I ran in pieces, because of patience, takes awhile!
st = 39
se = 39
for id in range(st,se+1):
#for id in range(1,len(deptht)):  # use this to run all at once
    print id
    for x in range(0,len(BX)):
        for y in range(0,len(BY)):
            if deptht[id] - depth[id,y,x] > 0.01:
                vosaline[0,id,y,x] = salinity_function(depth[id,y,x])
39
In [44]:
nemo.history = """
[2013-10-31 17:45:55] Created.
[2013-12-02 14:02:00] Changed to 428m max depth.
[2013-12-02 14:02:00] Updated to CF-1.6 conventions and netCDF4 with zlib compression.
[2013-12-31 12:15:00] Changed to smoothed JdF mouth.
[2014-01-24 16:49:00] Changed to variable depths.
"""
In [45]:
nc_tools.check_dataset_attrs(nemo)
In [46]:
nemo.close()

And do it all again for Temperature

In [47]:
nemo = NC.Dataset('SS2_SoG0318_1y_temperature_nomask.nc', 'w', zlib=True)
# dataset attributes
nc_tools.init_dataset_attrs(
    nemo, 
    title='Temperature Initial Condition', 
    notebook_name='PrepareTS', 
    nc_filepath='../../../nemo-forcing/initial_strat/SS2_SoG0318_1y_temperature_nomask.nc',
    comment='Temperature from S4-1 September 2003, STRATOGEM Cruise 03-18; used at all grid points')
# dimensions
nemo.createDimension('x',len(BX))
nemo.createDimension('y',len(BY))
nemo.createDimension('deptht',size = len(depth_out))
nemo.createDimension('time_counter', None)
# variables
nav_lat = nemo.createVariable('nav_lat', 'float32', ('y','x'))
nav_lat.long_name = 'Latitude'
nav_lat.units = 'degrees_north'
nav_lat = lat
nav_lon = nemo.createVariable('nav_lon', 'float32', ('y','x'))
nav_lon.long_name = 'Longitude'
nav_lon.units = 'degrees_east'
nav_lon = lon
deptht = nemo.createVariable('deptht', 'float32', ('deptht'))
deptht.long_name = 'Vertical T Levels'
deptht.units = 'm'
deptht.positive = 'down'
deptht.valid_range = np.array((4., 428.))
deptht = depth_out
time_counter = nemo.createVariable('time_counter', 'float32', ('time_counter'))
time_counter.units = 'seconds since 2003-09-10 12:27:00'
time_counter.long_name = 'Time axis'
time_counter.calendar = 'noleap' 
time_counter[0] = 1
votemper = nemo.createVariable('votemper', 'float32', 
                               ('time_counter','deptht','y','x'))
votemper.units = 'degC'
votemper.long_name = 'Temperature'  
file format: NETCDF4
Conventions: CF-1.6
title: Temperature Initial Condition
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://bitbucket.org/salishsea/tools/src/tip/PrepareTS.ipynb
references: https://bitbucket.org/salishsea/nemo-forcing/src/tip/initial_strat/SS2_SoG0318_1y_temperature_nomask.nc
history: [2014-01-25 09:45:39] Created netCDF4 zlib=True dataset.
comment: Temperature from S4-1 September 2003, STRATOGEM Cruise 03-18; used at all grid points

Assign temperature values:

In [48]:
votemper[0,0,:,:] = temp[0]
for id in range(1,len(deptht)):
    votemper[0,id,:,:] = temperature_function(deptht[id])
print votemper[0,:,15,20]
[ 14.04310036  14.04074955  14.00984955  13.94554806  13.80773544
  13.16073608  12.15436363  11.5372057   11.31503582  11.19921017
  11.08286476  10.98369312  10.9316988   10.90688896  10.86747456
  10.81823158  10.7899332   10.68832302  10.5691185   10.45423412
  10.31723785   9.91960526   9.94678497   9.67330265   9.44548798
   9.44641876   9.78447723   9.8582325    9.81627083   9.62160778
   9.54338932   9.41520977   9.41339302   9.30292892   9.38290977
   9.45765018   9.46421432   9.45713997   9.45899963   9.45899963]
In [50]:
# I ran in pieces, because of patience, takes awhile!
st = 16
se = 39
for id in range(st,se+1):
#for id in range(1,len(deptht)):   # use this to do it all at once
    print id
    for x in range(0,len(BX)):
        for y in range(0,len(BY)):
            if deptht[id] - depth[id,y,x] > 0.01:
                votemper[0,id,y,x] = temperature_function(depth[id,y,x])
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
In [51]:
nemo.history = """
[2013-10-31 17:45:55] Created.
[2013-12-02 14:04:00] Changed to 428m max depth.
[2013-12-02 14:04:00] Updated to CF-1.6 conventions and netCDF4 with zlib compression.
[2013-12-31 12:15:00] Changed to smoothed JdF mouth.
[2014-01-24 16:49:00] Changed to variable depths.
"""
In [52]:
nc_tools.check_dataset_attrs(nemo)
In [53]:
nemo.close()
In [ ]: