calc_tauw.m 4.61 KB
Newer Older
Pierre Cazenave's avatar
Pierre Cazenave committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
function [tauw] = calc_tauw(z0,ncfile,writenetcdf) 

% DESCRIPTION:
%    Function used to calculate Wave Bed Shear Stress, tauw, from the
%    bottom orbital velocity (Ubot) and the bottom wave period (TmBot). 
%    This function is currently set up to extract these values from SWAN
%    output. Currently temperature, salinity and z0 are set to be constant
%    accross the domain. Using actual data could improve accuaracy.
%    
%    REQUIRES FVCOM_TOOLBOX to be in your MATLAB path
% 
%
% INPUT 
%   z0           = bed roughness length (z0=d50/12) [mm]
%   <=={possibly use constant value 0.05 Soulsby pg. 49}==>
%   
%   ncfile       = netcdf file containing Ubot & TmBot
%                  that was created from swan2netcdf.m
%   
%   writenetcdf  = accepts either true or false 
%                  if 'True' write tauw and fw to netcdf
%   
%    
%   
%
% OUTPUT:
%    tauw = array containing values of tauw at all time for each node in 
%    model domain
%
% EXAMPLE USAGE
%    tauw = calc_tauw(z0,'skg4.3.nc',true); 
%    
%
% Author(s):  
%    Eric Holmes (University of Massachusetts Dartmouth)
%
% Revision history
%    Initially created 09-24-2010
%   
%==============================================================================

%------------------------------------------------------------------------------
% Convert z0 to meters
%------------------------------------------------------------------------------
z0 = z0/1000;

%------------------------------------------------------------------------------
% Set Ubot & TmBot from netCDF file
%------------------------------------------------------------------------------
nc = netcdf(ncfile,'w');

nctime=ncdim('time',nc);
nctime = size(nctime);
nctime=nctime(1,1);
ncnode=ncdim('node',nc);
node = size(ncnode);
node = node(1,1);
Ubot = ncvar('Ubot',nc);
Ubot = Ubot.*ones(nctime,node);
TmBot = ncvar('TmBot',nc);
TmBot = TmBot.*ones(nctime,node);



%------------------------------------------------------------------------------
% Set Generic Values to Salinity and Temperature
% This is temporary, it would be better to use actual salinity and
% temperature data.
%------------------------------------------------------------------------------
vectorsize=size(Ubot);
T=zeros(vectorsize);
S=zeros(vectorsize);

T(:,:)=15;
S(:,:)=30;

%------------------------------------------------------------------------------
% Call Kinematic Viscosity Routine and Density Routine
%------------------------------------------------------------------------------
nu = SW_Kviscosity(T,S); % 
rho = SW_Density(T,S); %


%------------------------------------------------------------------------------
% Calculate fwr & fws according to Soulsby pg. 77-79
%------------------------------------------------------------------------------

%----CONSTANTS-----------------------------------------------------------------
%ks = z0*30.;                    % Nikuradse roughness
A = Ubot.*TmBot/(2.*pi);        % semi-orbital excursion
% r = A/ks;                     % relative roughness
Rw = Ubot.*A./nu;                 % wave Reynolds
%------------------------------------------------------------------------------
fwr = 1.39*(A/z0).^(-0.52);     % case in which flow is rough turbulent 
                                % Soulsby (62a)

                                % Smooth Cases
for i = 1:vectorsize(1,1)
    for j = 1:vectorsize(1,2)
        if(Rw(i,j) > 5e5)
            B=0.0521; N=0.187;  % case in which flow is smooth turbulent
        else
            B=2; N=0.5;         % case in which flow is laminar
        end;
    end
end
                                
fws = B*Rw.^(-N);               % smooth bed friction factor Soulsby (63)

%------------------------------------------------------------------------------
% Choose wave friction factor for current conditions
%------------------------------------------------------------------------------
fw = zeros(vectorsize);

for i = 1:vectorsize(1,1)
    for j = 1:vectorsize(1,2)
        fw(i,j) = max(fwr(i,j),fws(i,j));   
                                % wave friction factor
    end;
end;

tauw = 0.5*rho.*fw.*Ubot.*Ubot;    % wave shear stress
%tauw(isnan(tauw)) = 0;

%------------------------------------------------------------------------------
% Write Values of tauw to ncfile
%------------------------------------------------------------------------------
if (writenetcdf == true)
        
    nc = netcdf(ncfile,'w');
    
    nc{'tauw'} = ncfloat('time','node');
    nc{'tauw'}.long_name = 'Bed Shear Stress';
    nc{'tauw'}.units     = 'N/m^2';
    
    for i=1:vectorsize(1,1)
        nc{'tauw'}(i,1:vectorsize(1,2)) = tauw(i,1:vectorsize(1,2));
    end
    
    close(nc);
    
end;









end