IRF Example Marine

From CSDMS
Revision as of 14:05, 26 October 2009 by Huttone (talk | contribs) (Before and after code to demonstrate IRF (marine example))
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

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()
        
    #----------------------------------------------------------------------