Commit 894fc35a 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 5ef1daf6 db5fb66b
......@@ -19,6 +19,8 @@ function example_FVCOM_river()
% 2016-02-18 TODO: This function uses add_var_FVCOM_river which in turn
% uses the old netCDF toolbox for MATLAB. That function needs to be
% updated to use the built in netCDF routines in MATLAB.
% 2017-03-31 add_var_FVCOM_river now uses the native MATLAB netCDF
% routines, so this script should work fine.
%
%==============================================================================
......@@ -49,22 +51,20 @@ RiverName = {'tstRiver'};
write_FVCOM_river(RiverFile,RiverName,time,flux,temp,salt,RiverInfo1,RiverInfo2)
% add sediment to the file
if exist('netcdf') == 2
VarName = 'fine_sand';
VarLongName = 'concentration of fine sand';
VarUnits = 'kgm^-3';
VarData = .333*sedload;
add_var_FVCOM_river(RiverFile,VarName,VarLongName,VarUnits,VarData)
VarName = 'fine_sand';
VarLongName = 'concentration of fine sand';
VarUnits = 'kgm^-3';
VarData = .333*sedload;
add_var_FVCOM_river(RiverFile,VarName,VarLongName,VarUnits,VarData)
VarName = 'coarse_silt';
VarLongName = 'concentration of coarse silt';
VarUnits = 'kgm^-3';
VarData = .333*sedload;
add_var_FVCOM_river(RiverFile,VarName,VarLongName,VarUnits,VarData)
VarName = 'coarse_silt';
VarLongName = 'concentration of coarse silt';
VarUnits = 'kgm^-3';
VarData = .333*sedload;
add_var_FVCOM_river(RiverFile,VarName,VarLongName,VarUnits,VarData)
VarName = 'fine_silt';
VarLongName = 'concentration of fine silt';
VarUnits = 'kgm^-3';
VarData = .333*sedload;
add_var_FVCOM_river(RiverFile,VarName,VarLongName,VarUnits,VarData)
end
VarName = 'fine_silt';
VarLongName = 'concentration of fine silt';
VarUnits = 'kgm^-3';
VarData = .333*sedload;
add_var_FVCOM_river(RiverFile,VarName,VarLongName,VarUnits,VarData)
......@@ -25,6 +25,7 @@ function write_FVCOM_bedflag(bedflag,filename,mytitle)
% 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.
% 2017-03-29 Write the flag as a float as FVCOM expects.
%
%==========================================================================
......@@ -70,7 +71,7 @@ netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'), 'title', mytitle)
node_dimid=netcdf.defDim(nc, 'node', numel(bedflag));
% define variables and attributes
node_varid=netcdf.defVar(nc, 'bedflag', 'NC_INT', node_dimid);
node_varid=netcdf.defVar(nc, 'bedflag', 'NC_FLOAT', node_dimid);
netcdf.putAtt(nc,node_varid, 'long_name', 'bed deposition flag');
netcdf.putAtt(nc,node_varid, 'units', '-');
......
......@@ -375,7 +375,7 @@ for ii = 1:numvars
fprintf('NEW DATA... ')
end
tmp_time = [];
for i = 1:nt;
for i = 1:nt
tmp_time = [tmp_time, sprintf('%-026s', datestr(datenum(out_date(i, :)), 'yyyy-mm-dd HH:MM:SS.FFF'))];
end
netcdf.putVar(ncout, varid, tmp_time)
......
......@@ -8,10 +8,14 @@ function write_FVCOM_sediment(sediment, filename)
% for FVCOM
%
% INPUT
% sediment = struct with nodal data in the following fields:
% sediment = struct with nodal data in the following mandatory fields:
% - tau_ce
% - tau_cd
% - tau_cd (this appears to be unused but required at launch)
% - erate
% and the following optional field(s):
% - <sediment_class>_bedfrac
% If the optional field(s) is included, the namelist option
% SEDIMENT_PARAMETER_TYPE = "non-uniform" can be used.
% filename = filename to which to write the outputs.
%
% OUTPUT:
......@@ -63,18 +67,30 @@ netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'), 'history', ...
node_dimid = netcdf.defDim(nc, 'node', nVerts);
% define variables and attributes
erate_varid = netcdf.defVar(nc, 'erate', 'NC_INT', node_dimid);
erate_varid = netcdf.defVar(nc, 'erate', 'NC_DOUBLE', 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);
tau_ce_varid = netcdf.defVar(nc, 'tau_ce', 'NC_DOUBLE', 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);
tau_cd_varid = netcdf.defVar(nc, 'tau_cd', 'NC_DOUBLE', 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');
extra_fields = setdiff(fieldnames(sediment), {'tau_ce', 'tau_cd', 'erate'});
for f = 1:length(extra_fields)
sediment_name = strrep(extra_fields{f}, '_bedfrac', '');
if ftbverbose
fprintf('Adding extra variable: %s\n', sediment_name)
end
eval(sprintf('%s_varid = netcdf.defVar(nc, ''%s'', ''NC_DOUBLE'', node_dimid);\n', sediment_name, sediment_name))
eval(sprintf('netcdf.putAtt(nc, %s_varid, ''long_name'', ''Fraction of %s in surface layer'');\n', sediment_name, sediment_name))
eval(sprintf('netcdf.putAtt(nc, %s_varid, ''units'', ''-'')\n', sediment_name))
eval(sprintf('netcdf.putAtt(nc, %s_varid, ''grid'', ''fvcom_grid'')\n', sediment_name))
end
% end definitions
netcdf.endDef(nc);
......
......@@ -27,8 +27,11 @@ function [dstar] = ST_Dstar(d,varargin)
%
%==============================================================================
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
fprintf('\nbegin : %s\n', subname)
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
% constants
grav = 9.8106; %g
......@@ -71,4 +74,6 @@ else
dstar = ([grav*(s-1)./(nu.^2)]).^(1/3).*d;
end
fprintf('end : %s\n', subname)
if ftbverbose
fprintf('end : %s\n', subname)
end
function [erate] = ST_erate(d,varargin)
function erate = ST_erate(d,varargin)
% Calculate erosion rate in kg/(m^2-s)
%
% function [erate] = ST_erate(d,varargin)
......@@ -34,11 +34,13 @@ function [erate] = ST_erate(d,varargin)
%
%==============================================================================
subname = 'ST_erate';
%fprintf('\n')
%fprintf(['begin : ' subname '\n'])
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
% constants
% constants
grav = 9.8106; %g
T = 10; %T (C)
S = 35; %S (PSU)
......@@ -70,6 +72,6 @@ wset = ST_wset(d,'temperature',T,'salinity',S,'sdens',sdens);
% calculate erosion rate
erate = 2.666e-4*wset*1000. - 2.51e-9*sdens;
%fprintf(['end : ' subname '\n'])
if ftbverbose
fprintf('end : %s\n', subname)
end
\ No newline at end of file
......@@ -16,7 +16,11 @@ function ST_example
% Revision history
%
%==============================================================================
close all;
global ftbverbose
ftbverbose = false;
close all
fprintf(' phi class d(mm) Dstar wset(mm/s) taucr (Pa) erate x1e-3(kg/(m^2-s))\n')
i = 0;
for phi=-8:11
......
function [d] = ST_phi2d(phi)
function d = ST_phi2d(phi)
% Convert sediment grain size from phi to d (m)
%
% function [d] = ST_phi2d(phi)
......@@ -22,9 +22,11 @@ function [d] = ST_phi2d(phi)
%
%==============================================================================
%subname = 'ST_phid2';
%fprintf('\n')
%fprintf(['begin : ' subname '\n'])
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
%------------------------------------------------------------------------------
% calculate d and convert to m
......@@ -33,4 +35,6 @@ function [d] = ST_phi2d(phi)
d = 2^(-phi);
d = d*.001;
%fprintf(['end : ' subname '\n'])
if ftbverbose
fprintf('end : %s\n', subname)
end
\ No newline at end of file
......@@ -29,8 +29,11 @@ function [taucr] = ST_taucr(d,varargin)
% 2017-03-27 Add support for matrices.
%==============================================================================
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
fprintf('\nbegin : %s\n', subname)
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
% constants
grav = 9.8106; %g
......@@ -74,4 +77,6 @@ else
taucr = theta_cr*grav*(sdens-dens)*d;
end
fprintf('end : %s\n', subname)
if ftbverbose
fprintf('end : %s\n', subname)
end
function [Sclass] = wentworth(phi)
function Sclass = ST_wentworth(phi)
% Report wentworth class of a particular grain size phi
%
% function wentworth(phi)
......@@ -22,6 +22,12 @@ function [Sclass] = wentworth(phi)
%
%==============================================================================
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
ClassNames = {'boulder','cobble','pebble','granule','very coarse sand', ...
'coarse sand','medium sand','fine sand','very fine sand','coarse silt', ...
'medium silt','fine silt','very fine silt','coarse clay','medium clay','fine clay'};
......@@ -30,3 +36,7 @@ pts = find(phi-ClassLbound > 0);
ClassIndex = pts(end);
Sclass = char(ClassNames{ClassIndex});
%fprintf('class of phi = %f is: %s\n',phi,class);
if ftbverbose
fprintf('end : %s\n', subname)
end
\ No newline at end of file
......@@ -28,11 +28,13 @@ function [wset] = ST_wset(d,varargin)
%
%==============================================================================
subname = 'ST_wset';
%fprintf('\n')
%fprintf(['begin : ' subname '\n'])
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
% constants
% constants
grav = 9.8106; %g
T = 10; %T (C)
S = 35; %S (PSU)
......@@ -70,7 +72,12 @@ dens = SW_Density(T,S);
dstar = ST_Dstar(d,'temp',T,'sal',S,'sdens',sdens);
% calculate wset
wset = (nu/d)*( sqrt(10.36^2 + 1.049*(dstar^3)) - 10.36);
if ismatrix(d)
wset = (nu./d).*( sqrt(10.36^2 + 1.049*(dstar.^3)) - 10.36);
else
wset = (nu/d)*( sqrt(10.36^2 + 1.049*(dstar^3)) - 10.36);
end
%fprintf(['end : ' subname '\n'])
if ftbverbose
fprintf('end : %s\n', subname)
end
\ No newline at end of file
......@@ -19,6 +19,7 @@ function rho = SW_Density(T,S)
%
% AUTHOR:
% Mostafa H. Sharqawy 12-18-2009, MIT (mhamed@mit.edu)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
......@@ -36,8 +37,17 @@ function rho = SW_Density(T,S)
% substance, 1996.
% UPDATED 09-23-2010 modified to now handle matrices and commented out
% range checking.
%
% Revision history
% 2017-03-31 Tidy up the code a bit.
%=========================================================================
global ftbverbose
[~, subname] = fileparts(mfilename('fullpath'));
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
......@@ -50,12 +60,16 @@ end
[mt,nt] = size(T);
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
if (ms~=mt) || (ns~=nt)
error('check_stp: S & T must have same dimensions')
end
% CHECK THAT S & T ARE WITHIN THE FUNCTION RANGE
vectorsize=size(S);
s = nan(vectorsize);
rho_w = nan(vectorsize);
D_rho = nan(vectorsize);
rho = nan(vectorsize);
for i = 1:vectorsize(1,1)
for j = 1:vectorsize(1,2)
% if T(i,j)>180 | T(i,j)<0
......@@ -75,5 +89,8 @@ for i = 1:vectorsize(1,1)
D_rho(i,j) = b1*s(i,j) + b2*s(i,j)*T(i,j) + b3*s(i,j)*T(i,j)^2 + b4*s(i,j)*T(i,j)^3 + b5*s(i,j)^2*T(i,j)^2;
rho(i,j) = rho_w(i,j) + D_rho(i,j);
end
end;
end
if ftbverbose
fprintf('end : %s\n', subname)
end
\ No newline at end of file
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