Demonstration using the NetCDF4-Python library to compute bottom stress and bottom velocity (1 meter above bottom) from a triangular grid ocean model (FVCOM) via OPeNDAP. The results are stored in a new NetCDF4 file.
NECOFS (Northeastern Coastal Ocean Forecast System) is run by groups at the University of Massachusetts Dartmouth and the Woods Hole Oceanographic Institution, led by Drs. C. Chen, R. C. Beardsley, G. Cowles and B. Rothschild. Funding is provided to run the model by the NOAA-led Integrated Ocean Observing System and the State of Massachusetts.
NECOFS is a coupled numerical model that uses nested weather models, a coastal ocean circulation model, and a wave model. The ocean model is a volume-mesh model with horizontal resolution that is finer in complicated regions. It is layered (not depth-averaged) and includes the effects of tides, winds, and varying water densities caused by temperature and salinity changes.
In [ ]:
from pylab import *
import netCDF4
In [ ]:
def z0tocd(z0=3.3546e-04,zr=1.0):
""" Compute the drag coefficient CD based on roughness height z0 and distance above bottom zr"""
kappa = 0.4
cd=(kappa*ones_like(zr)/log(zr/z0))**2
return cd
In [ ]:
def cdtoz0(cd=2.5e-3,zr=1.0):
""" Compute the roughness height z0 based on drag coefficient CD and distance above bottom zr"""
kappa = 0.4
z0 = zr/(exp(kappa*ones_like(cd)/sqrt(cd)))
return z0
In [ ]:
def w100(w=0.1+0j,z0=3.35e-04,zr=1):
""" Compute the velocity 1 meter above bottom and friction velocity
from velocity measured at height zr above bottom.
Keyword arguments
-----------------
w : east velocity component + j*north velocity component (m/s) [complex]
z0 : roughness height = kb/30 (m)
zr : height above bottom for input velocity "w" (m)
Returns
-------
ustar : friction velocity (m/s) [complex]
w : velocity 1 mab (m/s) [complex]
"""
cd = z0tocd(z0,zr)
ustar = sqrt(cd)*w
kappa = 0.4
ur = abs(ustar)/kappa*log(ones_like(zr)/z0)
wbot = w*ur/(abs(w)+1e-16)
return ustar,wbot
In [ ]:
# Input FVCOM Dataset: DAP Data URL
url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc'
# Open DAP
nci = netCDF4.Dataset(url)
nciv =nci.variables
nciv.keys()
In [ ]:
# read connectivity array
nv = nciv['nv'][:].T - 1
In [ ]:
# print info on velocity varibble
print nciv['u']
nt,nz,nele = shape(nciv['u'])
node = len(nciv['h'])
In [ ]:
# OUTPUT: create NetCDF4 file with deflation on variables
url_out = '/usgs/data2/rsignell/models/fvcom/massbay_forecast_wbot.nc'
nco = netCDF4.Dataset(url_out, 'w',clobber=True)
ncov = nco.variables
# chunking and deflation parameters
chunk_lon=512
chunk_time=1
sigdigits=6
# create dimensions
nco.createDimension('nele', nele)
nco.createDimension('node', node)
nco.createDimension('three', 3)
nco.createDimension('time', None)
# create variables
timeo = nco.createVariable('time', 'f4', ('time'))
ho = nco.createVariable('h', 'f4', ('node'))
nvo = nco.createVariable('nv', 'i4', ('three', 'nele'))
lonco = nco.createVariable('lonc', 'f4', ( 'nele'))
latco = nco.createVariable('latc', 'f4', ( 'nele'))
lono = nco.createVariable('lon', 'f4', ( 'node'))
lato = nco.createVariable('lat', 'f4', ( 'node'))
ubot = nco.createVariable('ubot', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
vbot = nco.createVariable('vbot', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
bustr = nco.createVariable('bustr', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
bvstr = nco.createVariable('bvstr', 'f4', ('time', 'nele'),
zlib=True,least_significant_digit=sigdigits,chunksizes=[chunk_time,chunk_lon])
# write variable attributes
timeo.units=nciv['time'].units
ho.units=nciv['h'].units
lono.units=nciv['lon'].units
lato.units=nciv['lat'].units
lonco.units=nciv['lonc'].units
latco.units=nciv['latc'].units
ubot.units=nciv['u'].units
vbot.units=nciv['v'].units
bustr.units='N/m2'
bvstr.units='N/m2'
# write data with no time dimension
lonco[:]=nciv['lonc'][:]
latco[:]=nciv['latc'][:]
lono[:]=nciv['lon'][:]
lato[:]=nciv['lat'][:]
nvo[:]=nciv['nv'][:]
ho[:]=nciv['h'][:]
In [ ]:
# specify bottom layer, but handle case where there is just 1 layer in input file
if shape(nciv['siglay'])[0]==1:
ilayer = 0
else:
ilayer = -1
In [ ]:
# neither z0 or cd is saved in this FVCOM output, so just use canonical bottom roughness kb=0.5 cm
kb=0.005
z0=kb/30.
rho = 1025. # density plays a small role in stress, so just specify as constant here
In [ ]:
# loop through time, writing each 2D or 3D field to output file
for itime in range(nt):
print '%4.1f percent done' % (float(itime)/nt*100.)
# get current at layer [0 = surface, -1 = bottom]
zr = (nciv['siglay'][ilayer,:]-nciv['siglev'][ilayer,:])*(nciv['h'][:]+nciv['zeta'][itime,:])
u = nciv['u'][itime, ilayer, :]
v = nciv['v'][itime, ilayer, :]
# average nodes to get bottom layer thicknesses at faces (velocity points)
zr_face = mean(zr[nv],axis=1)
# create complex velocity from components
w=u+1j*v
# compute bottom friction velocity and velocity at 1 mab
ustar,wbot = w100(w,z0,zr_face)
# compute bottom stress from friction velocity
bstr = rho * ustar * abs(ustar)
# write bottom velocity and stress components to output file
ubot[itime,:]=wbot.real
vbot[itime,:]=wbot.imag
bustr[itime,:]=bstr.real
bvstr[itime,:]=bstr.imag
In [ ]:
nci.close() # close input file (not stringly necessary for OPeNDAP, but good practice)
nco.close() # close output file
In [ ]: