Estimate time step at each node
[Mobj] = function estimate_ts(Mobj)
DESCRIPTION:
Calculate barotropic time step
INPUT
Mobj = matlab mesh object
OUTPUT:
Mobj = matlab structure containing mesh time step
EXAMPLE USAGE
Mobj = estimate_ts(Mobj)
Author(s):
Geoff Cowles (University of Massachusetts Dartmouth)
Revision history
2012-07-14 Add great circle approximation if only provided with
latitude and longitudes. Also add arguments to the function to define
current velocity and tidal amplitudes.
==============================================================================0001 function [Mobj] = estimate_ts(Mobj,u,zeta) 0002 0003 % Estimate time step at each node 0004 % 0005 % [Mobj] = function estimate_ts(Mobj) 0006 % 0007 % DESCRIPTION: 0008 % Calculate barotropic time step 0009 % 0010 % INPUT 0011 % Mobj = matlab mesh object 0012 % 0013 % OUTPUT: 0014 % Mobj = matlab structure containing mesh time step 0015 % 0016 % EXAMPLE USAGE 0017 % Mobj = estimate_ts(Mobj) 0018 % 0019 % Author(s): 0020 % Geoff Cowles (University of Massachusetts Dartmouth) 0021 % 0022 % Revision history 0023 % 2012-07-14 Add great circle approximation if only provided with 0024 % latitude and longitudes. Also add arguments to the function to define 0025 % current velocity and tidal amplitudes. 0026 % 0027 %============================================================================== 0028 0029 subname = 'estimate_ts'; 0030 global ftbverbose 0031 if(ftbverbose) 0032 fprintf('\n') 0033 fprintf(['begin : ' subname '\n']) 0034 end; 0035 0036 %------------------------------------------------------------------------------ 0037 % Set constants 0038 %------------------------------------------------------------------------------ 0039 0040 g = 9.81; %gravitational acceleration 0041 %u = 3.0; %u-velocity 0042 %zeta = 11.0; %tide amp 0043 0044 if(~Mobj.have_bath) 0045 error('can''t estimate the time step without bathymetry') 0046 end; 0047 0048 %------------------------------------------------------------------------------ 0049 % Compute the time step estimate 0050 %------------------------------------------------------------------------------ 0051 if Mobj.have_xy 0052 x = Mobj.x; 0053 y = Mobj.y; 0054 else 0055 % Will convert to metres when calculating element edge length 0056 x = Mobj.lon; 0057 y = Mobj.lat; 0058 end 0059 h = max(Mobj.h,0); 0060 tri = Mobj.tri; 0061 nVerts = Mobj.nVerts; 0062 nElems = Mobj.nElems; 0063 0064 ts = ones(nVerts,1)*1e9; 0065 lside = zeros(nVerts,1); 0066 for i=1:nElems 0067 n1 = tri(i,1); 0068 n2 = tri(i,2); 0069 n3 = tri(i,3); 0070 nds = [n1 n2 n3]; 0071 % Check whether we have x and y values and use great circle 0072 % approximations if we don't. 0073 if Mobj.have_xy 0074 lside(i) = sqrt( (x(n1)-x(n2))^2 + (y(n1)-y(n2))^2); 0075 else 0076 lside(i) = haversine(x(n1),y(n1),x(n2),y(n2)); 0077 end 0078 dpth = max(h(nds))+zeta; 0079 dpth = max(dpth,1); 0080 ts(nds) = min(ts(nds),lside(i)/(sqrt(g*dpth) + u)); 0081 end; 0082 if(ftbverbose); fprintf('minimum time step: %f seconds\n',min(ts)); end; 0083 Mobj.ts = ts; 0084 Mobj.have_ts = true; 0085 0086 if(ftbverbose) 0087 fprintf(['end : ' subname '\n']) 0088 end 0089 0090 function [km]=haversine(lat1,lon1,lat2,lon2) 0091 % Haversine function to calculate first order distance measurement. Assumes 0092 % spherical Earth surface. Lifted from: 0093 % 0094 % http://www.mathworks.com/matlabcentral/fileexchange/27785 0095 % 0096 R = 6371000; % Earth's mean radius in metres 0097 delta_lat = lat2 - lat1; % difference in latitude 0098 delta_lon = lon2 - lon1; % difference in longitude 0099 a = sin(delta_lat/2)^2 + cos(lat1) * cos(lat2) * ... 0100 sin(delta_lon/2)^2; 0101 c = 2 * atan2(sqrt(a), sqrt(1-a)); 0102 km = R * c; % distance in metres 0103