Commit 5c2da156 authored by Pierre Cazenave's avatar Pierre Cazenave

Merge branch 'swaldman3-rebased_rory' into dev

parents 1f066668 ceb638b0
function M = change_shallow_bathy(M, min_depth)
% Deepens shallow nodes by setting a minimum depth and making nodes that
% are too shallow the mean of the surrouding deep nodes. Loop until all
% nodes are deeper than the minuimum depth required.
%
% function [M]=change_shallow_bathy(M, min_depth)
%
% DESCRIPTION:
%
% INPUT:
% M = Mesh object
% min_depth = the minimum depth of nodes
%
% OUTPUT:
% Nested Mesh object with altered bathymetry.
%
% EXAMPLE USAGE:
%
%
% Author(s):
% Rory O'Hara Murray (Marine Scotland Science)
%
% Revision history:
% 2014 sometime - first version
%
%==========================================================================
count = 0;
h=-99; % just to start
while sum(h==-99)>0
count=count+1;
disp(['iteration # ' num2str(count)])
I = M.h<=min_depth; % find shallow nodes as depth is +ve down.
If = find(I);
h = [];
% loop through all shallow nodes with depth<=min_depth
for ii=1:length(If)
% find elements surrounding the shallow node
test = [];
for jj=1:3
test = [test; find(M.tri(:,jj)==If(ii))];
end
% find nodes for all these elements surrounding the shallow node
% (the surrounding nodes)
nodes = unique(M.tri(test,:));
htmp = M.h(nodes);
bla = htmp>min_depth; % find the nodes that are deeper than min_depth
if sum(bla)>0
h(ii) = mean(htmp(bla));% make the shallow node the mean of the surrounding nodes deeper than min_depth
else
h(ii) = -99; % id no deep surroundign nodes then make it -99 and try again.
end
end
M.h(I) = h; % save changes to the Mobj array.
end
end
\ No newline at end of file
......@@ -51,6 +51,12 @@ function [Mobj] = read_sms_mesh(varargin)
% does this anyway.
% 2016-07-28 Fix behaviour if grid has no open boundaries so we can rely
% on have_strings existing in either case.
% [the next few revisions are listed out of order because of rebasing
% a branch that had been separate for a long time]
% 2014-05-29 Changed the way the header is read and skipped (ROM).
% 2014-05-29 Changed the way the nodestrings are read, taking into
% account the possibility that SMS adds exatra 'name' number to each
% nodestring after the -ve indicator (ROM).
%
%==============================================================================
......@@ -183,11 +189,12 @@ lat = zeros(nVerts,1);
ts = zeros(nVerts,1);
% Skip the header
% C = textscan(fid, '%s', nHeader + 1,'Headerlines',nHeader); % why does it skip 3 lines when only two are in the header? Didn't use to throw errors though...
% if Meshname is more than one word the above would fail... Using
% headerlines as below is more safe proof.
for ii=1:nHeader
lin = fgetl(fid);
end
% Read the triangulation table
C = textscan(fid, '%s %d %d %d %d %d', nElems,'Headerlines',nHeader);
C = textscan(fid, '%s %d %d %d %d %d', nElems);
tri(:, 1) = C{3};
tri(:, 2) = C{4};
tri(:, 3) = C{5};
......@@ -215,29 +222,18 @@ end
% Build array of all the nodes in the nodestrings.
C = textscan(fid, '%s %d %d %d %d %d %d %d %d %d %d', nStrings);
allNodes = cell2mat(C(2:end))';
allNodes(allNodes == 0) = [];
% Add a new field to Mobj with all the nodestrings as a cell array.
nodeStrings = find(allNodes < 0);
read_obc_nodes = cell(1, length(nodeStrings));
for nString = 1:sum(allNodes(:) < 0)
if nString == 1
read_obc_nodes{nString} = abs(allNodes(1:nodeStrings(nString)));
else
read_obc_nodes{nString} = ...
abs(allNodes(nodeStrings(nString - 1) + 1: ...
nodeStrings(nString)));
end
% Check for closed nodestrings (which we don't really want).
if read_obc_nodes{nString}(1) == read_obc_nodes{nString}(end)
% Drop the end one.
warning('Closed node string found. Opening it.')
read_obc_nodes{nString} = read_obc_nodes{nString}(1:end - 1);
end
end
startp1 = 10*ceil(nodeStrings./10)+1;
ns_range = [[1; startp1(1:end-1)], nodeStrings];
if nStrings > 0
have_strings = true;
% Add a new field to Mobj with all the nodestrings as a cell array.
read_obc_nodes = cell(1, length(nodeStrings));
for nString = 1:size(ns_range,1)
read_obc_nodes{nString} = abs(allNodes(ns_range(nString,1):ns_range(nString,2)));
end
end
have_lonlat = false;
......@@ -376,6 +372,10 @@ assert(isfield(Mobj, 'x'), 'No coordinate data provided. Check your inputs and t
% Make a depth array for the element centres.
Mobj.hc = nodes2elems(h, Mobj);
% Add element spherical coordinates too.
Mobj.lonc = nodes2elems(lon, Mobj);
Mobj.latc = nodes2elems(lat, Mobj);
%--------------------------------------------------------------------------
% Add the Coriolis values
%--------------------------------------------------------------------------
......
......@@ -6,7 +6,7 @@
%
% Rory O'Hara Murray, 2013-07-01
%
function [Mobj TMD_ConList] = set_elevtide_tmd(Mobj, dates_MJD)
function [Mobj TMD_ConList] = set_elevtide_tmd(Mobj, dates_MJD, path_to_tmd)
%MyTitle = 'Julian FVCOM time series for open boundary from TPXO model using TMD';
......@@ -36,10 +36,7 @@ time = dates(1):1/24/6:dates(2);
time_MJD = time - datenum('1858-11-17 00:00:00');
model_file = 'DATA\Model_ES2008';
current_dir = pwd;
%tmp = which('tmd_tide_pred_2.m');
%a = strfind(tmp, 'tmd_tide_pred_2.m');
%cd(tmp(1:a-1))
cd([getenv('Hydro') '\Software\Matlab\TMD2.03\'])
cd(path_to_tmd);
[eta, TMD_ConList] = tmd_tide_pred_2(model_file, time, lat, lon, 'z');
cd(current_dir);
......
function strout = mjul2str(MJD,noyear)
% Convert a modified Julian day to a Matlab datestr style string
MJD = double(MJD);
mjul2matlab = 678942; %difference between modified Julian day 0 and Matlab day 0
if(~exist('noyear'))
strout = datestr(MJD+mjul2matlab);
......
......@@ -18,9 +18,38 @@
% [optional] axi = the axis - specify axis range
% [optional] pll = the axis
% [optional] grd = add gridlines - specify colour
% [optional] tit = add title - specify text
% [optional] leg = add legend - specify text as a cell array
% [optional] qui = add quiver vectors - specify a structure with the quiver information in (see examples)
%
%
% EXAMPLE USAGE
% plot_fvcom_field(Mobj, Mobj.zeta, 'fid', 1, 'cli', [0 100], 'gif', 'animation.gif', 'axi', [60000 70000 40000 50000])
% plot_fvcom_field(Mobj, Mobj.zeta, 'fid', 1, 'cli', [0 100], 'gif', 'animation.gif', 'axi', [60000 70000 40000 50000])
%
% Quiver vecotor example 1 (plot every other depth average velocity vector on unstructured grid):
% Q.X = PFOW.lonc(1:15:end);
% Q.Y = PFOW.latc(1:15:end);
% Q.U = PFOW.ua(1:15:end,:);
% Q.V = PFOW.va(1:15:end,:);
% plot_fvcom_field(PFOW, PFOW.ua(:,1:13), 'pll', 'qui', Q)
%
% Quiver vecotor example 2 (include vecotrs on an interpolated regular grid):
% Q.x = -4:0.01:-2;
% Q.y = 58:0.01:59;
%
% [Q.X1, Q.Y1] = meshgrid(Q.x, Q.y);
% Q.X = Q.X1(:); Q.Y = Q.Y1(:);
%
% Only use data from within the region of interpolation
% I = PFOW.lonc>Q.x(1) & PFOW.lonc<Q.x(end) & PFOW.latc>Q.y(1) & PFOW.latc<Q.y(end);
%
% for tt=1:13
% Fx = scatteredInterpolant(double(PFOW.lonc(I)), double(PFOW.latc(I)), double(PFOW.ua(I,tt)));
% Fy = scatteredInterpolant(double(PFOW.lonc(I)), double(PFOW.latc(I)), double(PFOW.va(I,tt)));
% Q.U(:,tt) = Fx(Q.X, Q.Y);
% Q.V(:,tt) = Fy(Q.X, Q.Y);
% end
% plot_fvcom_field(PFOW, PFOW.ua(:,1:13), 'pll', 'qui', Q)
%
% Author(s)
% Rory O'Hara Murray (Marine Scotland Science)
......@@ -60,6 +89,7 @@ axis_flag = false;
title_flag = false;
legend_text_flag = false;
quiver_flag = false;
quiver2_flag = false;
for ii=1:1:length(varargin)
keyword = lower(varargin{ii});
......@@ -92,6 +122,9 @@ for ii=1:1:length(varargin)
case 'qui'
quiver_flag = true;
quiverData = varargin{ii+1};
% case 'qu2'
% quiver2_flag = true;
% quiverData = varargin{ii+1};
end
end
......@@ -125,9 +158,17 @@ elseif size(plot_field,1)==size(x,1) % plot on nodes
end
end
if not(fig_flag)
if fig_flag
if fig.Type(1)=='f'
the_axes = axes;
elseif fig.Type(1)=='a'
the_axes = fig;
end
else
fig = figure;
the_axes = axes;
end
axes(the_axes);
for ii=1:size(plot_field,2)
if ishandle(fig)==0 break; end
......@@ -142,8 +183,12 @@ for ii=1:size(plot_field,2)
title(['time = ' datestr(double(M.time(ii))+MJD_datenum, 'HH:MM dd/mm/yyyy')])
end
if quiver_flag
% hold on
% quiver(quiverData.X, quiverData.Y, quiverData.U(:,:,ii), quiverData.V(:,:,ii), 'w');
% hold off
% elseif quiver2_flag
hold on
quiver(quiverData.X, quiverData.Y, quiverData.U(:,:,ii), quiverData.V(:,:,ii), 'k');
quiver(quiverData.X, quiverData.Y, quiverData.U(:,ii), quiverData.V(:,ii), 0.4, 'color', 0.99*[1 1 1])
hold off
end
......
......@@ -2,17 +2,22 @@
%
% Example Usage:
%
% sigma_frac = turbine_area_sigma(H, Ht, r, sigLay, plot_fig)
% sigma_frac = turbine_area_sigma(H, Ht, r, sigLay, plot_fig, subplot_info)
%
% Input Parameters: H - mean sea level (m)
% Ht - height of turbine hub above seabed (m)
% r - turbine rotor radius (m)
% sigLay - number of sigma layers (not levels) in the model
% plot_fig - flag to plot a figure (optional)
% plot_fig - optional; flag to plot a figure
% subplot_info - optional; if present should be an
% array containing the three parameters to subplot
% that should be used to put the figure into a
% subplot.
%
% Rory O'Hara Murray, 19-Nov-2014
% Simon Waldman, 2016.
%
function sigma_frac = turbine_area_sigma(H, Ht, r, sigLay, plot_fig)
function sigma_frac = turbine_area_sigma(H, Ht, r, sigLay, plot_fig, subplot_info)
assert(nargin >= 4, 'Not enough arguments.');
assert(isnumeric(sigLay) && sigLay - fix(sigLay) < eps, 'sigLay (4th parameter) must be an integer number of sigma layers.');
......@@ -20,17 +25,18 @@ assert(isnumeric(sigLay) && sigLay - fix(sigLay) < eps, 'sigLay (4th parameter)
if nargin<5
plot_fig = false;
end
if nargin<6
splot = false;
else
splot = true;
end
% Turbine and sigma layer parameters
elev = 0; % water elevation above/below MSL - change this to see how the sigma layer occupation fraction changes with the tide
depth = H + elev;
dT = depth - Ht; % turbine hub depth
dT = H - Ht; % turbine hub depth
assert(dT>r, 'Turbine will stick out of water');
dLay = depth./sigLay;
zLev = [0:-dLay:-depth]';
dLay = H./sigLay;
zLev = [0:-dLay:-H]';
% what sigma layer is the hub in?
drsl = zLev+dT; % depth of hub relative to each sigma level
......@@ -38,11 +44,15 @@ hub_sigma = sum(drsl>=0);
%% draw sigma levels / layers
if plot_fig
figure
if splot
subplot( subplot_info(1), subplot_info(2), subplot_info(3) )
else
figure
end
plot([-r r], zLev*[1 1])
xlabel('Distance (m)')
ylabel('Depth (m)')
title([num2str(depth, '%2.0f') ' m water depth'])
title([num2str(H, '%2.0f') ' m water depth'])
% draw rotor area
a=0;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment