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
==============================================================================0001 function [tauw] = calc_tauw(z0,ncfile,writenetcdf) 0002 0003 % DESCRIPTION: 0004 % Function used to calculate Wave Bed Shear Stress, tauw, from the 0005 % bottom orbital velocity (Ubot) and the bottom wave period (TmBot). 0006 % This function is currently set up to extract these values from SWAN 0007 % output. Currently temperature, salinity and z0 are set to be constant 0008 % accross the domain. Using actual data could improve accuaracy. 0009 % 0010 % REQUIRES FVCOM_TOOLBOX to be in your MATLAB path 0011 % 0012 % 0013 % INPUT 0014 % z0 = bed roughness length (z0=d50/12) [mm] 0015 % <=={possibly use constant value 0.05 Soulsby pg. 49}==> 0016 % 0017 % ncfile = netcdf file containing Ubot & TmBot 0018 % that was created from swan2netcdf.m 0019 % 0020 % writenetcdf = accepts either true or false 0021 % if 'True' write tauw and fw to netcdf 0022 % 0023 % 0024 % 0025 % 0026 % OUTPUT: 0027 % tauw = array containing values of tauw at all time for each node in 0028 % model domain 0029 % 0030 % EXAMPLE USAGE 0031 % tauw = calc_tauw(z0,'skg4.3.nc',true); 0032 % 0033 % 0034 % Author(s): 0035 % Eric Holmes (University of Massachusetts Dartmouth) 0036 % 0037 % Revision history 0038 % Initially created 09-24-2010 0039 % 0040 %============================================================================== 0041 0042 %------------------------------------------------------------------------------ 0043 % Convert z0 to meters 0044 %------------------------------------------------------------------------------ 0045 z0 = z0/1000; 0046 0047 %------------------------------------------------------------------------------ 0048 % Set Ubot & TmBot from netCDF file 0049 %------------------------------------------------------------------------------ 0050 nc = netcdf(ncfile,'w'); 0051 0052 Ubot = nc{'Ubot'}(:,:); 0053 TmBot = nc{'TmBot'}(:,:); 0054 0055 0056 0057 %------------------------------------------------------------------------------ 0058 % Set Generic Values to Salinity and Temperature 0059 % This is temporary, it would be better to use actual salinity and 0060 % temperature data. 0061 %------------------------------------------------------------------------------ 0062 vectorsize=size(Ubot); 0063 T=zeros(vectorsize); 0064 S=zeros(vectorsize); 0065 0066 T(:,:)=15; 0067 S(:,:)=30; 0068 0069 %------------------------------------------------------------------------------ 0070 % Call Kinematic Viscosity Routine and Density Routine 0071 %------------------------------------------------------------------------------ 0072 nu = SW_Kviscosity(T,S); % 0073 rho = SW_Density(T,S); % 0074 0075 0076 %------------------------------------------------------------------------------ 0077 % Calculate fwr & fws according to Soulsby pg. 77-79 0078 %------------------------------------------------------------------------------ 0079 0080 %----CONSTANTS----------------------------------------------------------------- 0081 %ks = z0*30.; % Nikuradse roughness 0082 A = Ubot.*TmBot/(2.*pi); % semi-orbital excursion 0083 % r = A/ks; % relative roughness 0084 Rw = Ubot.*A./nu; % wave Reynolds 0085 %------------------------------------------------------------------------------ 0086 fwr = 1.39*(A/z0).^(-0.52); % case in which flow is rough turbulent 0087 % Soulsby (62a) 0088 0089 % Smooth Cases 0090 for i = 1:vectorsize(1,1) 0091 for j = 1:vectorsize(1,2) 0092 if(Rw(i,j) > 5e5) 0093 B=0.0521; N=0.187; % case in which flow is smooth turbulent 0094 else 0095 B=2; N=0.5; % case in which flow is laminar 0096 end; 0097 end 0098 end 0099 0100 fws = B*Rw.^(-N); % smooth bed friction factor Soulsby (63) 0101 0102 %------------------------------------------------------------------------------ 0103 % Choose wave friction factor for current conditions 0104 %------------------------------------------------------------------------------ 0105 fw = zeros(vectorsize); 0106 0107 for i = 1:vectorsize(1,1) 0108 for j = 1:vectorsize(1,2) 0109 fw(i,j) = max(fwr(i,j),fws(i,j)); 0110 % wave friction factor 0111 end; 0112 end; 0113 0114 tauw = 0.5*rho.*fw.*Ubot.*Ubot; % wave shear stress 0115 %tauw(isnan(tauw)) = 0; 0116 0117 %------------------------------------------------------------------------------ 0118 % Write Values of tauw to ncfile 0119 %------------------------------------------------------------------------------ 0120 if (writenetcdf == true) 0121 0122 nc = netcdf(ncfile,'w'); 0123 0124 nc{'tauw'} = ncfloat('time','node'); 0125 nc{'tauw'}.long_name = 'Bed Shear Stress'; 0126 nc{'tauw'}.units = 'N/m^2'; 0127 0128 for i=1:vectorsize(1,1) 0129 nc{'tauw'}(i,1:vectorsize(1,2)) = tauw(i,1:vectorsize(1,2)); 0130 end 0131 0132 close(nc); 0133 0134 end; 0135 0136 0137 0138 0139 0140 0141 0142 0143 0144 end