Commit 3eb3adc2 authored by Ricardo Torres's avatar Ricardo Torres 💬
Browse files

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

parents 151330ef 2c8ae50b
function [Mobj] = add_stations_list(Mobj,Positions,Names,Dist)
% Add a set of stations at which FVCOM will output time series.
%
% [Mobj] = add_stations_list(Mobj,Positions,Names,Dist)
......
......@@ -427,9 +427,7 @@ for v = 1:length(fields)
% Now we've interpolated in space, we can interpolate the z-values
% to the sigma depths.
% Preallocate the output arrays
fvtempz = nan(nf, fz);
for pp = 1:nf
% Get the FVCOM depths at this node
tfz = sigma(oNodes(pp), :);
......
function write_FVCOM_bedflag(bedflag,filename,mytitle)
% Dump spatially-variable flag (bedflag) to FVCOM forcing file
function write_FVCOM_bedflag(bedflag,filename,mytitle)
% Write spatially-variable flag (bedflag) to FVCOM forcing file
%
% function write_FVCOM_bedflag(bedflag,filename,mytitle)
%
% DESCRIPTION:
% Generate a NetCDF file containing spatially variable bedflag for FVCOM
% Generate a netCDF file containing spatially variable bedflag for FVCOM
%
% INPUT
% bedflag = user defined bed flag (=0, no erosion/bedosition, =1, erosion/bedosition)
% INPUT
% bedflag = user defined bed flag (=0, no erosion/bedosition, =1, erosion/bedosition)
% on the nodes
% filename = filename to dump to
% mytitle = title of the case (set as global attribute)
% mytitle = title of the case (set as global attribute)
%
% OUTPUT:
% NetCDF file: filename
% netCDF file: filename
%
% EXAMPLE USAGE
% write_FVCOM_bedflag(bedflag, 'tst_bedflag.nc', 'no bedosition/erosion in Skagit river')
%
% Author(s):
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
% 2016-02-18 Updated the code to use the MATLAB netCDF routines.
%
%==============================================================================
% 2017-03-23 Add the supplied title to the generated netCDF file.
%
%==========================================================================
global ftbverbose
subname = 'write_FVCOM_bedflag';
......@@ -34,16 +34,16 @@ if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
%------------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Parse input arguments
%------------------------------------------------------------------------------
if(~exist('bedflag'))
%--------------------------------------------------------------------------
if ~exist('bedflag', 'var')
error('incorrect usage of gen_bedflag_file, must provide bedflag field')
end
if(~exist('filename'))
if ~exist('filename', 'var')
error('incorrect usage of gen_bedflag_file, must provide filename')
end
if(~exist('title'))
if ~exist('mytitle', 'var')
error('incorrect usage of gen_bedflag_file, must provide title field')
end
......@@ -53,16 +53,19 @@ if(nVerts == 0)
error('dimension of bedflag is 0, something is wrong ')
end;
%------------------------------------------------------------------------------
% Dump to bedflag NetCDF file
%------------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Dump to bedflag netCDF file
%--------------------------------------------------------------------------
if ftbverbose
fprintf('Dumping to bedflag NetCDF file: %s\m', filename);
fprintf('Dumping to bedflag netCDF file: %s\n', filename);
fprintf('Size of bedflag array: %d\n', nVerts);
end
nc = netcdf.create(filename, 'clobber');
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'), 'title', mytitle)
% define dimensions
node_dimid=netcdf.defDim(nc, 'node', numel(bedflag));
......@@ -81,5 +84,5 @@ netcdf.putVar(nc, node_varid, bedflag);
netcdf.close(nc)
if ftbverbose
fprintf('end : %s', subname)
fprintf('end : %s\n', subname)
end
function write_FVCOM_sediment(sediment, filename)
% Write spatially-varying sediment threshold data.
%
% function write_FVCOM_sediment(sediment, filename)
%
% DESCRIPTION:
% Generate a netCDF file containing spatially variable sediment propoerties
% for FVCOM
%
% INPUT
% sediment = struct with nodal data in the following fields:
% - tau_ce
% - tau_cd
% - erate
% filename = filename to which to write the outputs.
%
% OUTPUT:
% netCDF file: filename of the output netCDF file to generate.
%
% EXAMPLE USAGE
% write_FVCOM_sediment(sediment, 'tst_sediment.nc')
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
% 2017-03-27 New function to create a netCDF for the SEDIMENT_PARAMETER_FILE
% section in the model namelist.
%
%==========================================================================
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
%--------------------------------------------------------------------------
% Check inputs.
%--------------------------------------------------------------------------
assert(exist('sediment', 'var') == 1, 'incorrect usage of %s, must provide sediment struct', subname)
assert(exist('filename', 'var') == 1, 'incorrect usage of %s, must provide ''filename''', subname)
% Check dimensions
nVerts = numel(sediment.tau_cd);
assert(nVerts ~= 0, 'dimension of sediment data is 0, something is wrong')
%--------------------------------------------------------------------------
% Dump to netCDF file
%--------------------------------------------------------------------------
if ftbverbose
fprintf('Writing sediment data to netCDF file: %s\n', filename);
fprintf('nodes dimension: %d\n', nVerts);
end
nc = netcdf.create(filename, 'clobber');
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'), 'history', ...
sprintf('File created using %s from the MATLAB fvcom-toolbox', subname))
% define dimensions
node_dimid = netcdf.defDim(nc, 'node', nVerts);
% define variables and attributes
erate_varid = netcdf.defVar(nc, 'erate', 'NC_INT', node_dimid);
netcdf.putAtt(nc, erate_varid, 'long_name', 'bed erodibility constant');
netcdf.putAtt(nc, erate_varid, 'units', '-');
tau_ce_varid = netcdf.defVar(nc, 'tau_ce', 'NC_INT', node_dimid);
netcdf.putAtt(nc, tau_ce_varid, 'long_name', 'critical shear stress for erosion');
netcdf.putAtt(nc, tau_ce_varid, 'units', 'N/m^2');
tau_cd_varid = netcdf.defVar(nc, 'tau_cd', 'NC_INT', node_dimid);
netcdf.putAtt(nc, tau_cd_varid, 'long_name', 'critical shear stress for deposition');
netcdf.putAtt(nc, tau_cd_varid, 'units', 'N/m^2');
% end definitions
netcdf.endDef(nc);
% dump data
netcdf.putVar(nc, erate_varid, sediment.erate);
netcdf.putVar(nc, tau_ce_varid, sediment.tau_ce);
netcdf.putVar(nc, tau_cd_varid, sediment.tau_cd);
% close netCDF
netcdf.close(nc)
if ftbverbose
fprintf('end : %s\n', subname)
end
......@@ -4,7 +4,7 @@ function [dstar] = ST_Dstar(d,varargin)
% function [dstar] = ST_Dstar(d,varargin)
%
% DESCRIPTION:
% Convert grain size from d (m) to dimensionless D
% Convert grain size from d (m) to dimensionless D
%
% INPUT:
% d: sediment grain size in m
......@@ -16,20 +16,21 @@ function [dstar] = ST_Dstar(d,varargin)
% Dstar: nondimensional grain size
%
% EXAMPLE USAGE
% dstar = ST_Dstar(.0005,'temperature',10,'salinity',35,'sdens',2650)
% dstar = ST_Dstar(.0005,'temperature',10,'salinity',35,'sdens',2650)
%
% Author(s):
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
%
% 2017-03-27 Add support for matrices and some minor syntax changes.
%
%==============================================================================
subname = 'ST_Dstar';
%fprintf('\n')
%fprintf(['begin : ' subname '\n'])
[~, subname] = fileparts(mfilename('fullpath'));
fprintf('\nbegin : %s\n', subname)
% constants
% constants
grav = 9.8106; %g
T = 10; %T (C)
S = 35; %S (PSU)
......@@ -49,13 +50,12 @@ for i=1:2:length(varargin)-1
case 'sal'
S = varargin{i+1};
case 'sde'
sdens = varargin{i+1};
sdens = varargin{i+1};
otherwise
error(['Can''t understand value for:' keyword]);
end; %switch keyword
end;
% calculate nu
nu = SW_Kviscosity(T,S);
......@@ -63,7 +63,12 @@ nu = SW_Kviscosity(T,S);
dens = SW_Density(T,S);
% calculate dstar
s = sdens/dens;
dstar = ([grav*(s-1)/(nu^2)])^(1/3)*d;
if isscalar(dens)
s = sdens/dens;
dstar = ([grav*(s-1)/(nu^2)])^(1/3)*d;
else
s = sdens ./ dens;
dstar = ([grav*(s-1)./(nu.^2)]).^(1/3).*d;
end
%fprintf(['end : ' subname '\n'])
fprintf('end : %s\n', subname)
function [taucr] = ST_taucr(d,varargin)
% Calculate critical shear stress in Pascals
% Calculate critical shear stress in Pascals (equivalent to N/m^2)
%
% function [wset] = ST_taucr(d,varargin)
%
......@@ -13,26 +13,26 @@ function [taucr] = ST_taucr(d,varargin)
% [optional] 'sdens' = sediment density in kg/m^3 [default=2650]
%
% OUTPUT:
% taucr: critical shear stress in N/m^2
% taucr: critical shear stress in N/m^2
%
% EXAMPLE USAGE
% TCR = ST_taucr(.0005,'temperature',10,'salinity',35,'sdens',2650)
% TCR = ST_taucr(.0005,'temperature',10,'salinity',35,'sdens',2650)
%
% Author(s):
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% References
% Soulsby Dynamics of Marine Sands (SC77)
%
% Revision history
%
% 2017-03-27 Add support for matrices.
%==============================================================================
subname = 'ST_taucr';
%fprintf('\n')
%fprintf(['begin : ' subname '\n'])
[~, subname] = fileparts(mfilename('fullpath'));
fprintf('\nbegin : %s\n', subname)
% constants
% constants
grav = 9.8106; %g
T = 10; %T (C)
S = 35; %S (PSU)
......@@ -52,7 +52,7 @@ for i=1:2:length(varargin)-1
case 'sal'
S = varargin{i+1};
case 'sde'
sdens = varargin{i+1};
sdens = varargin{i+1};
otherwise
error(['Can''t understand value for:' keyword]);
end; %switch keyword
......@@ -63,14 +63,15 @@ end;
dens = SW_Density(T,S);
% calculate dstar
dstar = ST_Dstar(d,'temp',T,'sal',S,'sdens',sdens);
dstar = ST_Dstar(d, 'temp', T, 'sal', S, 'sdens', sdens);
% calculate theta_cr
theta_cr = (0.30/(1+1.2*dstar)) + 0.055*[1 - exp(-.020*dstar)];
% calculate theta_cr and then taucr
if ismatrix(dstar)
theta_cr = (0.30./(1+1.2.*dstar)) + 0.055.*[1 - exp(-.020.*dstar)];
taucr = theta_cr.*grav.*(sdens-dens).*d;
else
theta_cr = (0.30/(1+1.2*dstar)) + 0.055*[1 - exp(-.020*dstar)];
taucr = theta_cr*grav*(sdens-dens)*d;
end
% calculate taucr
taucr = theta_cr*grav*(sdens-dens)*d;
%fprintf(['end : ' subname '\n'])
fprintf('end : %s\n', subname)
function [array] = zero_to_nan(array)
function array = zero_to_nan(array)
%
% Replaces data outside error bands with an interpolated value
% Replaces values of zero with NaN.
%
[rr,cc]=size(array);
for ii=1:cc
ix = find((array(:,ii))==0.);
if ~isempty(ix)
for i = 1:length(ix);
array(ix(i),ii) = NaN;
end
end
end
% Author(s)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
% 2017-03-27 Removed the loops and used logical indexing instead.
array(array == 0) = nan;
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