Due to a shift in policy, from 0900 GMT on Wednesday 14th July 2021, we will be disabling ssh access to the server for external users. External users who wish to continue to access code repositories on the server will need to switch to using https. This can be accomplished in the following way: 1) On the repo on gitlab, use the clone dialogue and select ‘Clone with HTTPS’ to get the address of the repo; 2) From within the checkout of your repo run: $ git remote set-url origin HTTPS_ADDRESS. Here, replace HTTPS_ADDRESS with the address you have just copied from GitLab. Pulls and pushes will now require you to enter a username and password rather than using a ssh key. If you would prefer not to enter a password each time, you might consider caching your login credentials.

Commit ef72117c authored by Ricardo Torres's avatar Ricardo Torres 💬

Merge branch 'dev' of gitlab.ecosystem-modelling.pml.ac.uk:fvcom/fvcom-toolbox into dev

parents bdaa4913 d4e6c3d8
function [ maxCFLs, fighandle ] = show_max_CFL(grdFile, depFile, ncFile, extTS, coordtype, fig_flag)
%SHOW_MAX-CFL Function to find the max CFL encountered in each mesh element during an FVCOM model
%run. NB need enough RAM to load all u, v and zeta data from the given
%ncFile.
% Inputs: grdFile, depFile: casename_grd.dat and casename_dep.dat files
% from an FVCOM model
% ncFile : Output from FVCOM. Must include u, v, and zeta.
% extTS : External timestep at which to calculate CFL, in
% seconds.
% coordtype: 'lonlat' or 'cartesian'. If cartesian, coordinates
% should be in metres. If lonlat, they should be in
% degrees.
% fig_flag: Should a figure be plotted? true or false.
%
% The formulation of CFL used here is taken from the MIKE 3 manual
% Simon Waldman (Marine Scotland Science / Heriot-Watt University), March 2018
%% Check inputs
assert( nargin == 6, 'Wrong number of arguments.');
assert( exist(grdFile, 'file')==2, 'Cannot find file %s', grdFile );
assert( exist(depFile, 'file')==2, 'Cannot find file %s', depFile );
assert( exist(ncFile, 'file')==2, 'Cannot find file %s', ncFile );
assert( ischar(coordtype), 'Coord type must be ''lonlat'' or ''cartesian''.' );
assert( islogical( fig_flag ), 'fig_flag should be logical.' );
switch coordtype
case 'lonlat'
lonlat=true;
case 'latlon'
lonlat=true;
case 'cartesian'
lonlat=false;
otherwise
error('Coord type must be ''lonlat'' or ''cartesian''.');
end
% does the ncFile have velocities & water levels?
info = ncinfo(ncFile);
assert( any(strcmp('u', {info.Variables.Name})) && any(strcmp('v', {info.Variables.Name})) && ...
any(strcmp('zeta', {info.Variables.Name})), 'netCDF file must include the fields u, v and zeta.' );
%% Do the stuff.
disp('Reading mesh & bathymetry...');
M = read_fvcom_mesh( grdFile ); %NB this function puts lon and lat in the x and y fields.
M.h = read_fvcom_bath( depFile );
% calculate element depths from mean of nodes
M.hc = mean(M.h(M.tri),2);
disp('Reading U velocities...');
U = ncread(ncFile, 'u');
disp('Reading V velocities...');
V = ncread(ncFile, 'v');
disp('Reading surface elevations...');
Z = ncread(ncFile, 'zeta');
g = 9.81;
NumTS = size(U, 3);
NumEl = size(U, 1);
% Find water depths for each element at each TS. First calculate this for
% nodes, then average them for elements.
disp('Calculating depths...');
NodeDepths = repmat(M.h, 1, NumTS) + Z; %repmat because Z has a time dimension, M.h doesn't.
for n = 1:3
tmp(:,:,n) = NodeDepths(M.tri(:,n), :);
end
ElDepths = mean(tmp, 3); %this has dimensions of element x TS.
clear Z NodeDepths tmp; %save some memory
% Find minimum characteristic length for each element. This is approximated by the
% minimum edge length. It could be shorter with really long thin triangles,
% but the 30 degree internal angle minimum for FVCOM means we shouldn't have those.
disp('Calculating triangle sizes...');
CharLen = nan(NumEl, 1);
for e = 1:NumEl %for each element
xv = M.x(M.tri(e,:)); %vertices.
yv = M.y(M.tri(e,:));
%close the triangle by copying the first vertex to the end
xv(4) = xv(1);
yv(4) = yv(1);
%find the edge lengths
if lonlat
for a = 1:3
edges(a) = haversine(yv(a), xv(a), yv(a+1), xv(a+1));
end
else
edges = sqrt( diff(xv).^2 + diff(yv).^2 );
end
CharLen(e) = min(edges);
end
% For each element and TS, find the layer with the highest U and V magnitudes.
% Technically doing this with each component rather than the vector magnitude
% is wrong, but it'll usually be close and where wrong, it'll overestimate the CFL, so it's safe.
maxU = squeeze(max(abs(U), [], 2));
maxV = squeeze(max(abs(V), [], 2));
% Find CFL for each element at each TS
CFL = ( 2 .* sqrt( g .* ElDepths ) + maxU + maxV ) .* repmat( (extTS ./ CharLen), 1, NumTS );
%This is based on equation 6.1 on pg 33 of the MIKE hydrodynamic module
%manual (modified for using a single characteristic length rather than
%deltaX/deltaY)
% find the max over time for each element
MaxCFL = max(CFL, [], 2);
[val, I] = max(MaxCFL);
fprintf('Max CFL reached with an external timestep of %.2f secs is approx. %.3f, in Element %i.\n', extTS, val, I);
% find how long the timestep probably could go. We set the CFL to 0.8 and
% apply the same formula backwards.
TargetCFL = 0.8;
MaxTSs = repmat( (TargetCFL .* CharLen), 1, NumTS ) ./ ( 2 .* sqrt( g .* ElDepths ) + maxU + maxV );
%that's still per element per TS. We care about what the smallest is.
OverallMaxTS = min(min(MaxTSs));
fprintf('Max external timestep to reach CFL of %.1f with this mesh would be approx. %.2f seconds.\n', TargetCFL, OverallMaxTS );
% Optionally, plot figure of this
if fig_flag
CFLfig = figure;
p = patch();
p.Vertices = [M.x M.y];
p.Faces = M.tri;
p.CData = MaxCFL;
p.FaceColor = 'flat';
%p.EdgeColor = [0 0 0];
p.EdgeColor = 'none';
p.EdgeAlpha = 0.1;
p.LineWidth = 0.1;
cb = colorbar;
colormap(parula);
caxis([0 max(MaxCFL)]);
ylabel(cb, 'Max CFL encountered.');
axis equal;
% add to the plot markers for the worst n elements
n = 10; %number of cells to highlight
[ ~, WorstEls ] = sort(MaxCFL, 'descend');
WorstElsX = mean(M.x(M.tri(WorstEls(1:n),:)), 2); %finding element centroid coords.
WorstElsY = mean(M.y(M.tri(WorstEls(1:n),:)), 2);
hold on;
plot(WorstElsX, WorstElsY, 'or');
fprintf('Red circles on plot show the %i mesh elements with the highest CFL.\n', n);
end
%return values
maxCFLs = MaxCFL;
fighandle = CFLfig;
end
function [distm]=haversine(lat1,lon1,lat2,lon2)
% Haversine function to calculate first order distance measurement. Assumes
% spherical Earth surface. Lifted from:
%
% http://www.mathworks.com/matlabcentral/fileexchange/27785
% Could be done more accurately with Mapping Toolbox tools, but don't want
% to require that, and we don't need amazing accuracy.
lat1 = deg2rad(lat1);
lat2 = deg2rad(lat2);
lon1 = deg2rad(lon1);
lon2 = deg2rad(lon2);
R = 6371000; % Earth's mean radius in metres
delta_lat = lat2 - lat1; % difference in latitude
delta_lon = lon2 - lon1; % difference in longitude
a = sin(delta_lat/2)^2 + cos(lat1) * cos(lat2) * ...
sin(delta_lon/2)^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
distm = R * c; % distance in metres
end
\ No newline at end of file
......@@ -36,7 +36,7 @@ function [Mobj] = add_sponge_nodes_list(Mobj,SpongeList,SpongeName,SpongeRadius
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\m', subname)
fprintf('\nbegin : %s\n', subname)
end
% Do we want a figure showing how we're getting along?
......
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
......@@ -111,7 +111,12 @@ function Nested = find_nesting_region(conf, Mobj)
[~, subname] = fileparts(mfilename('fullpath'));
assert(isscalar(conf.power), 'conf.power should be scalar, applying to all open boundaries');
global ftbverbose
if isempty(ftbverbose)
ftbverbose=false;
end
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
......@@ -130,13 +135,13 @@ if any(conf.Nested_type == 3)
Nested.weight_node = cell(0);
end
if ftbverbose
figure(1)
clf
triplot(Nested.tri, Nested.x, Nested.y)
axis('equal', 'tight')
hold on
end
% if ftbverbose
% figure(1)
% clf
% triplot(Nested.tri, Nested.x, Nested.y)
% axis('equal', 'tight')
% hold on
% end
% Indices for the output cell arrays which are incremented for each nest
% level and each open boundary.
......@@ -238,9 +243,11 @@ for obc_idx = 1:Mobj.nObs
% setdiff
if ~isempty(truecandidate)
for dd=1:length(truecandidate)
Nested.read_obc_elems{end} = cat(2,Nested.read_obc_elems{end},truecandidate(dd));
Nested.weight_cell{end} = cat(2,Nested.weight_cell{end},Nested.weight_cell{end}(end));
end
Nested.read_obc_elems{end} = cat(2,Nested.read_obc_elems{end},truecandidate(dd));
if conf.Nested_type(obc_idx) ~= 1
Nested.weight_cell{end} = cat(2,Nested.weight_cell{end},Nested.weight_cell{end}(end));
end
end
end
if ftbverbose
fprintf('\n')
......
......@@ -58,11 +58,11 @@ while ~feof(fid)
if isempty(line) || strncmp(line, '!', 1) || ~ischar(line)
continue
end
% Clean up the input string to make matching a bit easier (trim
% whitespace and remove duplicate spaces in the keywords).
C = strtrim(regexpi(regexprep(line, '\s+', ' '), '=', 'split'));
switch lower(C{1})
case 'number of sigma levels'
nlev = str2double(C{2});
......
......@@ -26,6 +26,7 @@ function [Mobj] = read_sms_mesh(varargin)
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
% Rory O'Hara Murray (Marine Scotland Science)
% Simon Waldman (Marine Scotland Science / Heriot-Watt University)
%
% Revision history
%
......@@ -51,6 +52,18 @@ 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).
% 2018-05-16 Rewrote nodestring parsing. It's far less elegant, but now
% it still works if the number of nodes in a string is a multiple of 10.
% (SW)
% 2018-05-16 If we have bathymetry in the .2dm file *and* a separate
% bathymetry file provided, use the bathymetry in the file (with a
% warning) rather than ignoring it.
%
%==============================================================================
......@@ -183,11 +196,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};
......@@ -212,30 +226,51 @@ if max(isnan(tri(:))) == 1
error('%d NaNs in the h data', sum(isnan(tri(:))))
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);
%Read in nodestrings.
tmp = textscan(fid, ['%s' repmat('%d', 1, 12) '%*[^\n]'], nStrings, 'delimiter', ' ', 'MultipleDelimsAsOne', 1, 'CollectOutput', 1);
% this allows for up to 12 items on a NS line. It's normally 10, but can
% sometimes be 11. If we hit 12, something's changed in SMS's output.
% The second cell of the cell array returned by this should be a matrix of all the numeric
% values. Columns that don't have values in the file will contain 0.
mNSlines = tmp{2};
clear tmp;
% We'll work through the rows of this matrix and assemble the
% nodestring(s).
currentNSno = 1;
currentNS = [];
NSlengths = [];
for r = 1:size(mNSlines,1) %rows
for c = 1:size(mNSlines, 2) %columns
if mNSlines(r,c)==0 %we've run out of values on this line. Skip to next line.
break;
elseif mNSlines(r,c) < 0
%end of nodestring, marked by a negative node number.
%Append positive value to the NS, do end-of-NS stuff, then ignore
%rest of line.
currentNS = [currentNS abs(mNSlines(r,c))];
read_obc_nodes{currentNSno} = currentNS;
NSlengths = [NSlengths length(currentNS)];
currentNSno = currentNSno + 1;
currentNS = [];
break;
else
% append value to nodestring. If this is in a column higher
% than 10, raise a warning, as SMS doesn't usually do this.
currentNS = [currentNS mNSlines(r,c)];
if c > 10
warning('Longer lines than expected when parsing nodestrings; SMS output may not be as expected. Run with ftbverbose=true and check that the number and length of nodestrings found is as expected.');
end
end
end
end
if ftbverbose
a = sprintf('%d ', NSlengths);
fprintf('%i complete nodestrings found, of lengths %s. \n', currentNSno - 1, a);
clear a
end
if nStrings > 0
have_strings = true;
end
......@@ -268,38 +303,41 @@ fclose(fid);
bath_range = max(h) - min(h);
if have_bath
if bath_range == 0
fid = fopen(sms_bath, 'rt');
if fid < 0
error('file: %s does not exist', sms_bath);
else
if ftbverbose; fprintf('reading sms bathymetry from: %s\n', sms_bath); end
end
lin = fgetl(fid);
lin = fgetl(fid);
lin = fgetl(fid);
C = textscan(fid, '%s %d', 1);
nVerts_tmp = C{2};
C = textscan(fid, '%s %d', 1);
nElems_tmp = C{2};
if (nVerts - nVerts_tmp) * (nElems - nElems_tmp) ~= 0
fprintf('dimensions of bathymetry file do not match 2dm file\n')
fprintf('bathymetry nVerts: %d\n',nVerts_tmp)
fprintf('bathymetry nElems: %d\n',nElems_tmp)
error('stopping...')
end
lin = fgetl(fid);
lin = fgetl(fid);
lin = fgetl(fid);
lin = fgetl(fid); % extra one for the new format from SMS 10.1, I think
C2 = textscan(fid, '%f', nVerts);
h = C2{1};
have_bath = true;
clear C2
fclose(fid);
if bath_range ~= 0
warning(['Bathymetry is present in the .2dm file, but a bathymetry .dat file was also provided. '...
'The .dat file will be used, and the depth info in the .2dm will be ignored.']);
end
fid = fopen(sms_bath, 'rt');
if fid < 0
error('file: %s does not exist', sms_bath);
else
if ftbverbose; fprintf('reading sms bathymetry from: %s\n', sms_bath); end
end
lin = fgetl(fid);
lin = fgetl(fid);
lin = fgetl(fid);
C = textscan(fid, '%s %d', 1);
nVerts_tmp = C{2};
C = textscan(fid, '%s %d', 1);
nElems_tmp = C{2};
if (nVerts - nVerts_tmp) * (nElems - nElems_tmp) ~= 0
fprintf('dimensions of bathymetry file do not match 2dm file\n')
fprintf('bathymetry nVerts: %d\n',nVerts_tmp)
fprintf('bathymetry nElems: %d\n',nElems_tmp)
error('stopping...')
end
lin = fgetl(fid);
lin = fgetl(fid);
lin = fgetl(fid);
lin = fgetl(fid); % extra one for the new format from SMS 10.1, I think
C2 = textscan(fid, '%f', nVerts);
h = C2{1};
have_bath = true;
clear C2
fclose(fid);
elseif bath_range ~= 0
have_bath = true;
end
......@@ -342,8 +380,10 @@ if have_lonlat
Mobj.have_lonlat = have_lonlat;
Mobj.lon = lon;
Mobj.lat = lat;
Mobj.x = zeros(size(lon));
Mobj.y = zeros(size(lat));
if ~have_xy
Mobj.x = zeros(size(lon));
Mobj.y = zeros(size(lat));
end
% Add element spherical coordinates too.
Mobj.lonc = nodes2elems(lon, Mobj);
Mobj.latc = nodes2elems(lat, Mobj);
......@@ -352,8 +392,10 @@ if have_xy
Mobj.have_xy = have_xy;
Mobj.x = x;
Mobj.y = y;
Mobj.lon = zeros(size(x));
Mobj.lat = zeros(size(y));
if ~have_lonlat
Mobj.lon = zeros(size(x));
Mobj.lat = zeros(size(y));
end
% Add element cartesian coordinates too.
Mobj.xc = nodes2elems(x, Mobj);
Mobj.yc = nodes2elems(y, Mobj);
......@@ -376,6 +418,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);
......
......@@ -26,23 +26,26 @@ function write_FVCOM_river_nml(Mobj, nml_file, nc_file, vString)
% Z2=Mobj2.siglevz(Mobj2.river_nodes(r),2:end);
% RivDist(r,:)=(exp(Z1/5)-exp(Z2/5));
% end
%
% OUTPUT:
% Namelist for inclusion in the main FVCOM namelist (RIVER_INFO_FILE).
%
% EXAMPLE USAGE:
% write_FVCOM_river_nml(Mobj, 'casename_river.nml', 'casename_river.nc',RivDist)
% write_FVCOM_river_nml(Mobj, 'casename_river.nml', 'casename_river.nc', RivDist)
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
% Ricardo Torres (Plymouth Marine Laboratory)
%
% Revision history:
% 2013-03-21 - First version.
% 2013-08-16 - Add optional fourth argument of a string to supply as the
% RIVER_VERTICAL_DISTRIBUTION (e.g. 'uniform').
% 2013-10-16 - Fix the handling of the optional vString argument.
% 2018-03-28 - (RJT) Added the option of outputing float vertical distribution for river flows not just uniform for all rivers.
% 2018-03-28 - (RJT) Added the option of outputing float vertical
% distribution for river flows not just uniform for all rivers.
% 2018-05-08 - Fixed vertical distribution handling to support both the
% original functionality and the new float distribution.
%==========================================================================
subname = 'write_FVCOM_river_nml';
......@@ -57,17 +60,17 @@ nr = length(Mobj.river_nodes);
f = fopen(nml_file, 'w');
assert(f >= 0, 'Error writing to %s. Check permissions and try again.', nml_file)
for r = 1:nr
fprintf(f, ' &NML_RIVER\n');
fprintf(f, ' RIVER_NAME = ''%s'',\n', Mobj.river_names{r});
fprintf(f, ' RIVER_FILE = ''%s'',\n', nc_file);
fprintf(f, ' RIVER_GRID_LOCATION = %d,\n', Mobj.river_nodes(r));
% Build the vertical distribution string. Round to 15 decimal places so the
% unique check works (hopefully no one needs that many vertical layers...).
% Build the vertical distribution string. Round to 15 decimal places so
% the unique check works (hopefully no one needs that many vertical
% layers...).
vDist = roundn(abs(diff(Mobj.siglev)), -15);
if length(unique(vDist)) == 1
if length(unique(vDist)) == 1 || ~exist(vString, 'var')
vString = '''uniform''';
elseif nargin <= 3
vString = char();
......@@ -76,17 +79,16 @@ for r = 1:nr
for ii = 1:size(vDist)
vString = [vString, sprintf('%f ', vDist(ii))];
end
fprintf(f, ' RIVER_VERTICAL_DISTRIBUTION = %s\n', vString);
end
else % we have different vDist for each river
vString2=char();
for ii = 1:size(vString,2)
vString2 = [vString2, sprintf('%1.5f ', vString(r,ii))];
end
fprintf(f, ' RIVER_VERTICAL_DISTRIBUTION = %s\n', vString2);
else % we have different vDist for each river
vString = char();
for ii = 1:size(vString,2)
vString = [vString2, sprintf('%1.5f ', vString(r,ii))];
end
end
% Output whatever vertical distribution string we've got to the
% namelist file.
fprintf(f, ' RIVER_VERTICAL_DISTRIBUTION = %s\n', vString);
fprintf(f, ' /\n');
end
......
function RX1matrix=ComputeMatrixRx1_nodes(Z_w, Mobj)
% Computes rx1 matrix as in ComputeMatrixRx1_V2.m downloaded from
% Computes rx1 matrix as in ComputeMatrixRx1_V2.m downloaded from
% https://github.com/dcherian/tools
%
% function [RX1matrix] = ComputeMatrixRx1_nodes(Z_w,Mobj)
%
% function [RX1matrix] = ComputeMatrixRx1_nodes(Z_w, Mobj)
%
% DESCRIPTION:
% Calculates the hydrostatic consistency condition
% Calculates the hydrostatic consistency condition:
% r = abs(Sigma/H deltaxH/deltaSigma)
% this reflects the errors associated with horizontal pressure gradients calculations
% that are associated with steep bathyemtry and low vertical resolution.
% this reflects the errors associated with horizontal pressure gradients
% calculations that are associated with steep bathyemtry and low
% vertical resolution.
%
% INPUT
% Mobj = needs triangulation and mesh information
% table. read_sms_mesh provides everything it needs.
% INPUT
% Z_w = This is the sigmal layer vertical distribution in Z coordinates
%
% Mobj = needs triangulation and mesh information
% table. read_sms_mesh provides everything it needs.
%
%