Commit 4298320c authored by Rory O'Hara Murray's avatar Rory O'Hara Murray

Add the idealised tidal channel model by Rory O'Hara Murray

parent 7508f9e3
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
OBC Node Number = 145
1 8549 1
2 8550 1
3 8551 1
4 8552 1
5 8553 1
6 8554 1
7 8555 1
8 8556 1
9 8557 1
10 8605 1
11 8635 1
12 8661 1
13 8662 1
14 8663 1
15 8664 1
16 8665 1
17 8666 1
18 8667 1
19 8668 1
20 8669 1
21 8670 1
22 8647 1
23 8648 1
24 8619 1
25 8573 1
26 8512 1
27 8445 1
28 8373 1
29 8374 1
30 8284 1
31 8285 1
32 8286 1
33 8287 1
34 8288 1
35 8289 1
36 8290 1
37 8291 1
38 8292 1
39 8293 1
40 8294 1
41 8295 1
42 8296 1
43 8297 1
44 8298 1
45 8299 1
46 8300 1
47 8301 1
48 8375 1
49 8446 1
50 8513 1
51 8574 1
52 8620 1
53 8649 1
54 8650 1
55 8651 1
56 8671 1
57 8654 1
58 8655 1
59 8656 1
60 8657 1
61 8658 1
62 8659 1
63 8660 1
64 8632 1
65 8633 1
66 8589 1
67 8590 1
68 8591 1
69 8592 1
70 8593 1
71 8594 1
72 8595 1
73 8596 1
74 1 1
75 2 1
76 3 1
77 4 1
78 5 1
79 6 1
80 7 1
81 8 1
82 9 1
83 10 1
84 11 1
85 12 1
86 13 1
87 14 1
88 15 1
89 16 1
90 17 1
91 18 1
92 19 1
93 20 1
94 21 1
95 22 1
96 23 1
97 24 1
98 25 1
99 26 1
100 27 1
101 28 1
102 29 1
103 30 1
104 31 1
105 32 1
106 33 1
107 34 1
108 35 1
109 36 1
110 37 1
111 38 1
112 39 1
113 40 1
114 41 1
115 42 1
116 43 1
117 44 1
118 45 1
119 46 1
120 47 1
121 48 1
122 49 1
123 50 1
124 51 1
125 52 1
126 53 1
127 54 1
128 55 1
129 56 1
130 57 1
131 58 1
132 59 1
133 60 1
134 61 1
135 62 1
136 63 1
137 64 1
138 65 1
139 66 1
140 67 1
141 68 1
142 69 1
143 70 1
144 71 1
145 72 1
NUMBER OF SIGMA LEVELS = 11
SIGMA COORDINATE TYPE = UNIFORM
Sponge Node Number = 145
8549 150.000000 0.100000
8550 150.000000 0.100000
8551 150.000000 0.100000
8552 150.000000 0.100000
8553 150.000000 0.100000
8554 150.000000 0.100000
8555 150.000000 0.100000
8556 150.000000 0.100000
8557 150.000000 0.100000
8605 150.000000 0.100000
8635 150.000000 0.100000
8661 150.000000 0.100000
8662 150.000000 0.100000
8663 150.000000 0.100000
8664 150.000000 0.100000
8665 150.000000 0.100000
8666 150.000000 0.100000
8667 150.000000 0.100000
8668 150.000000 0.100000
8669 150.000000 0.100000
8670 150.000000 0.100000
8647 150.000000 0.100000
8648 150.000000 0.100000
8619 150.000000 0.100000
8573 150.000000 0.100000
8512 150.000000 0.100000
8445 150.000000 0.100000
8373 150.000000 0.100000
8374 150.000000 0.100000
8284 150.000000 0.100000
8285 150.000000 0.100000
8286 150.000000 0.100000
8287 150.000000 0.100000
8288 150.000000 0.100000
8289 150.000000 0.100000
8290 150.000000 0.100000
8291 150.000000 0.100000
8292 150.000000 0.100000
8293 150.000000 0.100000
8294 150.000000 0.100000
8295 150.000000 0.100000
8296 150.000000 0.100000
8297 150.000000 0.100000
8298 150.000000 0.100000
8299 150.000000 0.100000
8300 150.000000 0.100000
8301 150.000000 0.100000
8375 150.000000 0.100000
8446 150.000000 0.100000
8513 150.000000 0.100000
8574 150.000000 0.100000
8620 150.000000 0.100000
8649 150.000000 0.100000
8650 150.000000 0.100000
8651 150.000000 0.100000
8671 150.000000 0.100000
8654 150.000000 0.100000
8655 150.000000 0.100000
8656 150.000000 0.100000
8657 150.000000 0.100000
8658 150.000000 0.100000
8659 150.000000 0.100000
8660 150.000000 0.100000
8632 150.000000 0.100000
8633 150.000000 0.100000
8589 150.000000 0.100000
8590 150.000000 0.100000
8591 150.000000 0.100000
8592 150.000000 0.100000
8593 150.000000 0.100000
8594 150.000000 0.100000
8595 150.000000 0.100000
8596 150.000000 0.100000
1 150.000000 0.100000
2 150.000000 0.100000
3 150.000000 0.100000
4 150.000000 0.100000
5 150.000000 0.100000
6 150.000000 0.100000
7 150.000000 0.100000
8 150.000000 0.100000
9 150.000000 0.100000
10 150.000000 0.100000
11 150.000000 0.100000
12 150.000000 0.100000
13 150.000000 0.100000
14 150.000000 0.100000
15 150.000000 0.100000
16 150.000000 0.100000
17 150.000000 0.100000
18 150.000000 0.100000
19 150.000000 0.100000
20 150.000000 0.100000
21 150.000000 0.100000
22 150.000000 0.100000
23 150.000000 0.100000
24 150.000000 0.100000
25 150.000000 0.100000
26 150.000000 0.100000
27 150.000000 0.100000
28 150.000000 0.100000
29 150.000000 0.100000
30 150.000000 0.100000
31 150.000000 0.100000
32 150.000000 0.100000
33 150.000000 0.100000
34 150.000000 0.100000
35 150.000000 0.100000
36 150.000000 0.100000
37 150.000000 0.100000
38 150.000000 0.100000
39 150.000000 0.100000
40 150.000000 0.100000
41 150.000000 0.100000
42 150.000000 0.100000
43 150.000000 0.100000
44 150.000000 0.100000
45 150.000000 0.100000
46 150.000000 0.100000
47 150.000000 0.100000
48 150.000000 0.100000
49 150.000000 0.100000
50 150.000000 0.100000
51 150.000000 0.100000
52 150.000000 0.100000
53 150.000000 0.100000
54 150.000000 0.100000
55 150.000000 0.100000
56 150.000000 0.100000
57 150.000000 0.100000
58 150.000000 0.100000
59 150.000000 0.100000
60 150.000000 0.100000
61 150.000000 0.100000
62 150.000000 0.100000
63 150.000000 0.100000
64 150.000000 0.100000
65 150.000000 0.100000
66 150.000000 0.100000
67 150.000000 0.100000
68 150.000000 0.100000
69 150.000000 0.100000
70 150.000000 0.100000
71 150.000000 0.100000
72 150.000000 0.100000
This diff is collapsed.
function [PointID Distance] = fun_nearest2D(xloc, yloc, x, y)
%
% Find the nearest point in a 2D plane to a given location
%
% USAGE: [PointID Distance] = fun_nearest2D(xloc, yloc, x, y)
%
% xloc = x coordinate of location to find
% yloc = y coordinate of location to find
% x = x coordinates of all locations within 2D field
% y = y coordinates of all locations within 2D field
%
% Rory O'Hara Murray
%
for ii=1:length(xloc)
radvec = sqrt( (xloc(ii)-x).^2 + (yloc(ii)-y).^2);
[Distance(ii),PointID(ii)] = min(radvec);
end
return
\ No newline at end of file
% Plot an FVCOM field. This is somewhat similar to plot_field.m but is for
% postprocessing/viewing. It looks for the nv field included in fvcom
% output files. This function runs an animation if the field includes more
% than one time steps.
%
% plot_fvcom_field(Mobj, PlotField, 'fid', figure_id, 'cli', colour_lims, 'gif',
% filename, 'axi', axis_range, 'pll', 'grd', colour);
%
% INPUT
% Mobj = matlab mesh object with the following fields:
% - lon, lat, x, y = nodal positions (spherical and/or cartesian)
% - nv, tri = connectivity table (called either nv or tri)
% - [optional] time = Modified Julian Day time series
% PlotField = vertex-based field to plot
% [optional] fid = the fid of the figure to plot the field in - specify figure id
% [optional] cli = the colour limits to use - specify the limits
% [optional] gif = make an animated gif - specify filename
% [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])
%
% 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)
%
% Developments:
% 2014-05-22: Changed the way fig id is checked, not using 'exist' anymore.
% 2014-08-15: Added the axis command in
%
function [a] = plot_fvcom_field(M, plot_field, varargin)
MJD_datenum = datenum('1858-11-17 00:00:00');
% check to see if nv or tri should be used.
if isfield(M, 'nv')
nv = M.nv;
elseif isfield(M, 'tri')
nv = M.tri;
end
% check to see if a time variable is there or not
if isfield(M, 'time') %& size(M.time,1)>1
time_flag = true;
else
time_flag = false;
end
% defaults
clims = [min(plot_field(:)) max(plot_field(:))];
if clims(1)==clims(2)
clims(1)=clims(1)-0.1;
clims(2)=clims(2)+0.1;
end
gif = false;
grd = false;
plot_ll = false;
fig_flag = false;
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});
if length(keyword)~=3 || not(isstr(keyword))
continue
end
switch keyword
case 'fid' % id of a figure
fig = varargin{ii+1};
fig_flag = true;
case 'cli' % colour limits
clims = varargin{ii+1};
case 'gif' % make an animated gif
gif = true;
gif_filename = varargin{ii+1};
case 'axi' % axis
axis_flag = true;
axi = varargin{ii+1};
case 'grd' % grid lines
grd = true;
edgecolor = varargin{ii+1};
case 'pll'
plot_ll = true;
case 'tit'
title_flag = true;
fig_title = varargin{ii+1};
case 'leg'
legend_text_flag = true;
legend_text = varargin{ii+1};
case 'qui'
quiver_flag = true;
quiverData = varargin{ii+1};
if isfield(quiverData,'scale')==0 quiverData.scale = 1; end
if isfield(quiverData,'colour')==0 quiverData.colour =0.99*[1 1 1]; end
% case 'qu2'
% quiver2_flag = true;
% quiverData = varargin{ii+1};
end
end
if plot_ll
x = M.lon;
y = M.lat;
else
x = M.x;
y = M.y;
end
if not(axis_flag)
axi = [min(x) max(x) min(y) max(y)];
end
xE = x(nv)';
yE = y(nv)';
plot_field = squeeze(plot_field);
if size(plot_field,1)==size(nv,1) % plot on elements
if grd
patch_func = @(dummy) patch(xE, yE, dummy', 'edgecolor', edgecolor);
else
patch_func = @(dummy) patch(xE, yE, dummy', 'linestyle', 'none');
end
elseif size(plot_field,1)==size(x,1) % plot on nodes
if grd
patch_func = @(dummy) patch('Vertices',[x, y], 'Faces',nv, 'Cdata',dummy,'edgecolor', edgecolor,'facecolor','interp');
else
patch_func = @(dummy) patch('Vertices',[x, y], 'Faces',nv, 'Cdata',dummy,'linestyle','none','facecolor','interp');
end
end
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
a = patch_func(plot_field(:,ii));
c = colorbar;
if legend_text_flag set(get(c, 'ylabel'), 'string', legend_text); end
set(gca, 'clim', clims);
axis(axi)
if title_flag
title(fig_title)
elseif time_flag
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), quiverData.scale, 'color', quiverData.colour)
hold off
end
if gif
axis off
set(gcf, 'color', 'w')
frame = getframe(1);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if ii == 1;
imwrite(imind,cm,gif_filename,'gif', 'Loopcount',inf);
else
imwrite(imind,cm,gif_filename,'gif','WriteMode','append');
end
else
pause(0.01);
end
end
return
function M = read_fvcom_grid_data(filename, ll_only_flag)
if nargin<2
ll_only_flag = false;
end
M.lon = ncread(filename, 'lon');
M.lat = ncread(filename, 'lat');
M.nv = ncread(filename, 'nv');
M.h = ncread(filename, 'h');
M.lonc = ncread(filename, 'lonc');
M.latc = ncread(filename, 'latc');
if not(ll_only_flag)
M.x = ncread(filename, 'x');
M.y = ncread(filename, 'y');
M.xc = ncread(filename, 'xc');
M.yc = ncread(filename, 'yc');
end
M.lon(M.lon>180) = M.lon(M.lon>180) - 360;
M.lonc(M.lonc>180) = M.lonc(M.lonc>180) - 360;
M.nElems = length(M.lonc);
M.nNodes = length(M.latc);
return
\ No newline at end of file
% Load in idealised channel along channel velocity outputs and plot domain
% and a cross channel velocity section (animation through time), to
% visualise the output.
%
% Prerequisite functions should be bundled with this script
%
% Rory O'Hara Murray (r.murray@marlab.ac.uk)
% Marine Scotland Science
% 11 Dec 2019
%
clear, close all
filename1 = '../output/idealised_channel_02_0001.nc';
M1 = read_fvcom_grid_data(filename1);
M1.time = ncread(filename1, 'time');
M1.ua = ncread(filename1, 'ua', [1 1], [Inf length(M1.time)]);
M1.u = ncread(filename1, 'u', [1 1 1], [Inf Inf length(M1.time)]);
M1.zeta = ncread(filename1, 'zeta', [1 1], [Inf length(M1.time)]);
M1.siglay = ncread(filename1, 'siglay');
Ie1 = fun_nearest2D(5000,500,M1.xc, M1.yc);
In1 = fun_nearest2D(5000,500,M1.x, M1.y);
% define element centred depths
M1.hc = mean(M1.h(M1.nv),2);
for ii=1:size(M1.zeta,2)
M1.zetac(:,ii) = mean(M1.zeta(M1.nv,ii),2);
end
% find elements that lie within +/- 25 m of a cross channel transect line
IeT = find(M1.xc>=4975 & M1.xc<5025);
[B,I] = sort(M1.yc(IeT));
IeT = IeT(I);
% find nodes that lie within +/- 25 m of a along channel transect line
InTa = find(M1.y>=475 & M1.y<525 & M1.x>0 & M1.x<10000);
[B,I] = sort(M1.x(InTa));
InTa = InTa(I);
%% plot the model domain
plot_fvcom_field(M1, M1.h, 'tit', 'Idealised Channel domain and bathymetry')
axis equal
c = colorbar;
set(get(c, 'ylabel'), 'string', 'Depth (m)')
%% make an animation showing the velocity cross section
fig=figure('position', [48 70 1082 688]);
nSig = 10;
y = M1.yc(IeT) * ones(1,nSig);
for tt=length(M1.time)+[-floor(12.42*3):1:0]%200:300
if ishandle(fig)==0 break; end
u = squeeze(M1.u(IeT,:,tt));
%d = (ones(length(IeT),3)+M1.siglay(IeT,:)) .* (M1.hc(IeT)*ones(1,3));
d = -M1.siglay(IeT,:) .* ((M1.hc(IeT)+M1.zetac(IeT,tt))*ones(1,nSig));
subplot(4,1,[1 2]);
pcolor(y, d, u)
% axis([0 1000 0 40]), caxis([-2.5 2.5])
set(gca, 'ydir', 'reverse', 'ylim', [0 40], 'xlim', [0 1000], 'clim', [-1.5 1.5])
shading interp
c = colorbar;
set(get(c, 'ylabel'), 'string', 'Along channel velocity (m/s)')
xlabel('Distance across channel (m)')
ylabel('Depth (m)')
subplot(4,1,3)
plot(M1.x(InTa), M1.zeta(InTa,tt))
grid on
xlabel('Distance along channel (m)')
ylabel('Water elevation (m)')
set(gca, 'ylim', 1.3*[-1 1])
subplot(4,1,4);
plot(M1.time, M1.ua(Ie1,:), M1.time(tt), M1.ua(Ie1, tt), 'or')
grid on
set(gca, 'ylim', [-1.5 1.5])
xlabel('Time (days)')
ylabel('Depth mean velocity (m/s)')
pause(0.25)
end
% Create all idealised channel FVCOM model input files, including a spectral tidal boundary forcing file.
%
% Prerequisites: FVCOM matlab toolbox
%
% Rory O'Hara Murray (r.murray@marlab.ac.uk)
% Marine Scotland Science
% 11 Dec 2019
%
clear, close all
casename = 'idealised_channel_02';
inputdir = '../input/';
smsdir = inputdir;
global ftbverbose
ftbverbose = true; %print information to screen [true/false]
% read the Mesh from an SMS file
Mobj = read_sms_mesh('2dm',[smsdir casename '.2dm'],'project','false');