Commit 5175f32d authored by Pierre Cazenave's avatar Pierre Cazenave
Browse files

Update the read_sigma function to better deal with geometric distributions....

Update the read_sigma function to better deal with geometric distributions. Also clean up the code to read the parameters from the input file (this should be a bit more robust) and add a few checks the input is correct
parent a8f86d8b
......@@ -25,17 +25,21 @@ function Mobj = read_sigma(Mobj, sigmafile)
% Revision history
% 2013-01-08 Based on the code in show_sigma.m but instead of calculating
% sigma layers along a user-defined line, the depths are calculated for
% sigma layer along a user-defined line, the depths are calculated for
% each node in the unstructured grid.
% 2013-01-10 Added two new outputs to the returned Mobj (siglay and
% siglev). They're useful in write_FVCOM_tsobc.m.
% 2013-04-23 Add support for geometric sigma distributions as well as
% slightly more robust reading of and checks on the parameters in the
% input file. Also changed the way the uniform distribution is calculated
% (by using a P_SIGMA value of 1 and the sigma_geo.m function rather than
% fiddling around with ranges, although the output is the same).
subname = 'read_sigma';
global ftbverbose;
if ftbverbose
fprintf(['begin : ' subname '\n'])
fprintf(['\nbegin : ' subname '\n'])
fid = fopen(sigmafile,'r');
......@@ -48,22 +52,27 @@ while ~feof(fid)
if isempty(line) || strncmp(line, '!', 1) || ~ischar(line)
key = lower(line(1:3));
C = strtrim(regexpi(line, '=', 'split'));
switch key
case 'num'
% 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});
case 'sig'
case 'sigma coordinate type'
sigtype = C{2};
case 'du '
case 'sigma power'
sigpow = str2double(C{2});
case 'du'
du = str2double(C{2});
case 'dl '
case 'dl'
dl = str2double(C{2});
case 'min'
case 'min constant depth'
min_constant_depth = str2double(C{2});
case 'ku '
case 'ku'
ku = str2double(C{2});
case 'kl '
case 'kl'
kl = str2double(C{2});
case 'zku'
s = str2double(regexp(C{2}, ' ', 'split'));
......@@ -80,27 +89,51 @@ while ~feof(fid)
% Do some checks if we've got uniform or generalised coordinates to make
% sure the input is correct.
if strcmpi(sigtype, 'GENERALIZED') || strcmpi(sigtype, 'UNIFORM')
if numel(zku) ~= ku
warning('Number of zku values does not match the number specified in ku')
if numel(zkl) ~= kl
warning('Number of zkl values does not match the number specified in kl')
if ftbverbose
fprintf('nlev %d\n',nlev)
fprintf('sigtype %s\n',sigtype)
fprintf('du %d\n',du)
fprintf('dl %d\n',dl)
fprintf('min_constant_depth %f\n',min_constant_depth)
fprintf('ku %d\n',ku)
fprintf('kl %d\n',kl)
fprintf('zku %d\n',zku)
fprintf('zkl %d\n',zkl)
% Should be present in all sigma files.
fprintf('nlev\t%d\n', nlev)
fprintf('sigtype\t%s\n', sigtype)
% Only present in geometric sigma files.
if strcmpi(sigtype, 'GEOMETRIC')
fprintf('sigpow\t%d\n', sigpow)
% Only in the generalised or uniform sigma files.
if strcmpi(sigtype, 'GENERALIZED') || strcmpi(sigtype, 'UNIFORM')
fprintf('du\t%d\n', du)
fprintf('dl\t%d\n', dl)
fprintf('min_constant_depth\t%f\n', min_constant_depth)
fprintf('ku\t%d\n', ku)
fprintf('kl\t%d\n', kl)
fprintf('zku\t%d\n', zku)
fprintf('zkl\t%d\n', zkl)
% calculate the sigma distributions at each grid node
% Calculate the sigma distributions at each grid node
switch lower(sigtype)
case 'generalized'
z = sigma_gen(nlev, dl, du, kl, ku, zkl, zku, ...
Mobj.z(i), min_constant_depth);
case 'uniform'
z = 0:-1/double(nlev-1):-1;
z = sigma_geo(nlev, 1);
% z = 0:-1/double(nlev-1):-1;
case 'geometric'
z = sigma_geo(nlev, sigpow);
error('Can''t do that sigtype')
error('Don''t know sigtype %s (is it supported?)', sigtype)
% Create a siglay variable (i.e. midpoint in the sigma levels).
......@@ -115,4 +148,4 @@ Mobj.siglay = zlay;
if ftbverbose;
fprintf(['end : ' subname '\n'])
\ 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