IRF Example Marine
From CSDMS
The MATLAB code to calculate bottom stress:
% bottom_stress.m by Carl Friedrichs, 20 Oct 2009
% input: hourly wave orbital vel from wave_orb_vel.dat ascii vector
% output: bed stress taub.dat ascii vector
load ub.dat % hourly mean amplitude of orbital velocity (m/s)
fw=0.008; % wave friction factor (dimensionless)
rho=1028; % water density (kg/m^3)
taub=rho*fw/2*ub.^2; % bottom stress (Pascals)
figure(1); plot(ub); title('input: ub = hourly orbital velocity (m/s)'); % plot of input
figure(2); plot(taub); title('output: taub = hourly bottom stress (Pa)'); % plot of output
save taub.dat taub -ascii % one column ascii output file
The IRF Python code to calculate bottom stress:
# October 20, 2009
#
#-------------------------------------------------------------------
import numpy
#-------------------------------------------------------------------
class model_class():
def initialize(self):
self.fw = numpy.float32(0.008) # [meters]
self.rho = 1028.0 # [kg/m^3]
self.taub_list = []
self.DONE = False
#----------------------------------------------------------------------
def update(self, VERBOSE=False):
taub = self.rho * (self.fw/2) * self.ub**2
self.taub = taub
self.taub_list.append( taub )
if (VERBOSE):
print 'taub =', taub
#----------------------------------------------------------------------
def get_value(self, var_name):
if (var_name == 'taub'):
return self.taub
else:
print 'ERROR in get_value().'
#----------------------------------------------------------------------
def set_value(self, var_name, var):
if (var_name == 'ub'):
self.ub = var
else:
print 'ERROR in set_value().'
#----------------------------------------------------------------------
def finalize(self):
print 'Finished with simulation.'
#----------------------------------------------------------------------
def run(self):
self.initialize()
while not(self.DONE):
self.update()
self.finalize()
#----------------------------------------------------------------------
The MATLAB code to calculate orbital velocity:
% wave_orb_vel.m by Carl Friedrichs, 20 Oct 2009
% input: wave height in wave_height.dat ascii table (from Duck, NC)
% output: wave orbital velocity in ub.dat ascii vector
load wave_height.dat %yr,mo,dy,hr(EDT),ID,dir(deg),ht(m),period(s)
H=wave_height(:,7); % mean wave ht (m)
h=4; % water depth (m)
g=9.8; % accel of grav (m/s^2)
ub=0.5*(H/h)*sqrt(g*h); % hourly bottom orbital velocity (m/s)
figure(1); plot(H); title('input: H = hourly wave height (m)'); % plot of input
figure(2); plot(ub); title('output: ub = hourly orbital velocity (m/s)'); % plot of output
save ub.dat ub -ascii % one column ascii output file
The IRF Python code to calculate orbital velocity:
# October 20, 2009
#
#-------------------------------------------------------------------
import numpy
#-------------------------------------------------------------------
def test():
model = model_class()
model.run()
#-------------------------------------------------------------------
class model_class():
def initialize(self, file_name='/Users/peckhams/Desktop/Carl_Matlab/wave_height.dat' ):
self.file_unit = open( file_name, mode='r' )
self.h = numpy.float32(4) # [meters]
self.g = 9.8 # [m/s^2]
self.H = 0
self.ub = 0
self.DONE = False
self.ub_list = []
#----------------------------------------------------------------------
def update(self, VERBOSE=False):
line = self.file_unit.readline()
words = line.split()
if (len(words) >= 7):
g = self.g
h = self.h
self.H = float( words[6] )
self.ub = 0.5 * (self.H / self.h) * numpy.sqrt( g * h )
self.ub_list.append( self.ub )
if (VERBOSE):
print 'ub =', self.ub
self.DONE = (line == '')
#----------------------------------------------------------------------
def get_value(self, var_name):
if (var_name == 'H'):
return self.H
elif (var_name == 'ub'):
return self.ub
else:
print 'ERROR in get_value().'
#----------------------------------------------------------------------
def finalize(self):
self.file_unit.close()
print 'Finished with simulation.'
#----------------------------------------------------------------------
def run(self):
self.initialize()
while not(self.DONE):
self.update()
self.finalize()
#----------------------------------------------------------------------