Commit 466833ad authored by Pierre Cazenave's avatar Pierre Cazenave

Add support for the surface elevation open boundary node file _elj_obc.dat...

Add support for the surface elevation open boundary node file _elj_obc.dat (which is almost the same as tide_el.dat, but hey ho). Also make the error messages clearer and adjust the formatting statements to match what Dima has in his script
parent cedca07e
......@@ -42,28 +42,28 @@ if ftbverbose
fprintf('\n'); fprintf(['begin : ' subname '\n']);
end
%% _meanflow.dat
%% _meanflow.dat -- mean flow velocities at the open boundary elements (?).
% Create depth averaged velocity from the 3D velocity data in Mobj.
velocity = squeeze(mean(sqrt(Mobj.meanflow_u.^2 + Mobj.meanflow_v.^2), 2));
f = fopen([casename, '_meanflow.dat'], 'w');
if f < 0
error('Problem writing to .dat file. Check permissions and try again.')
error('Problem writing to _meanflow.dat file. Check permissions and try again.')
end
% Number of boundary nodes
fprintf(f, '%i\n', numel(Mobj.read_obc_nodes{1}));
fprintf(f, '%8d\n', Mobj.nObcNodes);
% Boundary node IDs
for i = 1:numel(Mobj.read_obc_nodes{1})
fprintf(f, '%i\n', Mobj.read_obc_nodes{1}(i));
for i = 1:Mobj.nObcNodes
fprintf(f, '%8d\n', Mobj.read_obc_nodes{1}(i));
end
% Sigma level distribution
s = '%i\t';
s = '%8d';
for ss = 1:length(Mobj.siglay)
if ss < length(Mobj.siglay)
s = [s, '%.4f\t'];
s = [s, '%8.4f'];
else
s = [s, '%.4f\n'];
s = [s, '%8.4f\n'];
end
end
for i = 1:numel(Mobj.read_obc_nodes{1})
......@@ -76,12 +76,12 @@ end
% Add the number of time steps
fprintf(f, '%i\n', nt);
s = '%.4f\n';
s = '%8.4f\n';
for ss = 1:nb
if ss < nb
s = [s, '%.4f\t'];
s = [s, '%8.4f'];
else
s = [s, '%.4f\n'];
s = [s, '%8.4f\n'];
end
end
for i = 1:size(velocity, 2)
......@@ -93,10 +93,10 @@ fclose(f);
%% _tide_node.dat -- nodes along the open boundaries.
f = fopen([casename, '_tide_node.dat'], 'w');
if f < 0
error('Problem writing to .dat file. Check permissions and try again.')
error('Problem writing to _tide_node.dat file. Check permissions and try again.')
end
% Boundary node IDs
fprintf(f, '%i\n', numel(Mobj.obc_nodes(Mobj.obc_nodes ~= 0)));
fprintf(f, '%8d\n', numel(Mobj.obc_nodes(Mobj.obc_nodes ~= 0)));
for j = 1:length(Mobj.read_obc_nodes); % number of boundaries
for i = 1:numel(Mobj.read_obc_nodes{j})
fprintf(f, '%8i\n', Mobj.read_obc_nodes{j}(i));
......@@ -108,14 +108,14 @@ fclose(f);
%% _tide_cell.dat -- elements which have two nodes on an open boundary.
f = fopen([casename, '_tide_cell.dat'], 'w');
if f < 0
error('Problem writing to .dat file. Check permissions and try again.')
error('Problem writing to _tide_cell.dat file. Check permissions and try again.')
end
if ~isfield(Mobj, 'read_obc_elements')
error('Missing list of boundary element IDs. Run find_boundary_elements and try again.')
end
% Boundary element IDs
ne = Mobj.nObcElements;
fprintf(f, '%i\n', ne);
fprintf(f, '%8d\n', ne);
for j = 1:length(Mobj.read_obc_nodes); % number of boundaries
for i = 1:numel(Mobj.read_obc_elements{j})
fprintf(f, '%8i\n', Mobj.read_obc_elements{j}(i));
......@@ -124,10 +124,10 @@ end
fclose(f);
%% _tide_el.dat
%% _tide_el.dat -- surface elevations with time.
f = fopen([casename, '_tide_el.dat'], 'w');
if f < 0
error('Problem writing to .dat file. Check permissions and try again.')
error('Problem writing to _tide_el.dat file. Check permissions and try again.')
end
% Boundary node IDs
if ~isfield(Mobj, 'surfaceElevation')
......@@ -139,12 +139,12 @@ end
[nb, nt] = size(Mobj.surfaceElevation);
s = '%i\t';
s = '%8d';
for ss = 1:nb
if ss < nb
s = [s, '%.4f\t'];
s = [s, '%8.4f'];
else
s = [s, '%.4f\n'];
s = [s, '%8.4f\n'];
end
end
......@@ -155,16 +155,15 @@ end
fclose(f);
%% _tide_uv.dat -- boundary velocities
% The format here is pretty funky. According to Dima's wrf_elj_obc.m
% The format here is pretty funky. According to Dima's wrt_elj_obc.m
% script, for each time step, there's lines of depth averaged u and v
% followed by all the u components at each vertical level, then all the v
% components at each vertical level. All lines are prefixed with the
% current time in seconds relative to the start of the model (or mean
% flow? God knows).
% followed by all the u and v components at each vertical level. All lines
% are prefixed with the current time in seconds relative to the start of
% the model (or mean flow time series? God knows).
f = fopen([casename, '_tide_uv.dat'], 'w');
if f < 0
error('Problem writing to .dat file. Check permissions and try again.')
error('Problem writing to _tide_uv.dat file. Check permissions and try again.')
end
% Number of elements in the boundaries.
......@@ -178,13 +177,13 @@ nz = length(Mobj.siglay);
% Create a format string for the each time step plus the number of boundary
% elements.
s = '%i\t';
s = '%8d';
for ss = 1:ne
if ss < ne
s = [s, '%.4f\t'];
s = [s, '%8.4f'];
else
s = [s, '%.4f\n'];
s = [s, '%8.4f\n'];
end
end
......@@ -209,8 +208,25 @@ end
fclose(f);
%% _elj_obc.dat -- surface elevation time series at open boundary nodes.
% This is almost identical to _tide_el.dat but lacks the time stamp in the
% first column.
f = fopen([casename, '_elj_obc.dat'], 'w');
if f < 0
error('Problem writing to _elj_obc.dat file. Check permissions and try again.')
end
nt = size(Mobj.surfaceElevation, 2);
for t = 1:nt
fprintf(f, '%8.4f', Mobj.surfaceElevation(:, t));
fprintf(f, '\n');
end
if ftbverbose
fprintf('end : %s \n', subname)
fprintf('end : %s\n', subname)
end
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