read_sigma.m 4.66 KB
Newer Older
1
function Mobj = read_sigma(Mobj, sigmafile)
2
% Read an FVCOM sigma layers file and output z values into Mobj.
3 4 5 6 7 8 9 10 11 12 13
% 
% Mobj = read_sigma(Mobj, sigmafile)
% 
% DESCRIPTION:
%   Read a sigma file and calculate the sigma layer depths
% 
% INPUT:
%   Mobj:       Mesh object which must contain Mobj.h (depths).
%   sigmafile : Full path to an FVCOM sigma.dat file.
% 
% OUTPUT:
14 15 16 17 18
%   Mobj:       Mesh object with four new fields:
%               - siglayz and siglevz: contain depths of the sigma layers
%               and levels at each grid node.
%               - siglay and siglev: the sigma layer and levels in the
%               range 0 to -1.
19 20
% 
% EXAMPLE USAGE:
21
%   Mobj = read_sigma(Mobj, 'sigma.dat')
22 23 24 25 26 27
% 
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
% 
% Revision history
%   2013-01-08 Based on the code in show_sigma.m but instead of calculating
Pierre Cazenave's avatar
Pierre Cazenave committed
28
%   sigma layers along a user-defined line, the depths are calculated for
29
%   each node in the unstructured grid.
30 31
%   2013-01-10 Added two new outputs to the returned Mobj (siglay and
%   siglev). They're useful in write_FVCOM_tsobc.m.
32 33 34 35 36
%   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).
37 38 39 40 41

subname = 'read_sigma';

global ftbverbose;
if ftbverbose
42
    fprintf(['\nbegin : ' subname '\n'])
43 44 45 46 47 48 49 50 51 52 53 54
end

fid = fopen(sigmafile,'r');
if(fid  < 0)
    error(['file: ' sigmafile ' does not exist']);
end

while ~feof(fid)
    line = fgetl(fid);
    if isempty(line) || strncmp(line, '!', 1) || ~ischar(line)
        continue
    end
55 56 57 58 59 60 61

    % 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'
62
            nlev = str2double(C{2});
63
        case 'sigma coordinate type'
64
            sigtype = C{2};
65 66 67
        case 'sigma power'
            sigpow = str2double(C{2});
        case 'du'
68
            du = str2double(C{2});
69
        case 'dl'
70
            dl = str2double(C{2});
71
        case 'min constant depth'
72
            min_constant_depth = str2double(C{2});
73
        case 'ku'
74
            ku = str2double(C{2});
75
        case 'kl'
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
            kl = str2double(C{2});
        case 'zku'
            s = str2double(regexp(C{2}, ' ', 'split'));
            zku = zeros(ku, 1);
            for i = 1:ku
                zku(i) = s(i);
            end
        case 'zkl'
            s = str2double(regexp(C{2}, ' ', 'split'));
            zkl = zeros(kl, 1);
            for i = 1:kl
                zkl(i) = s(i);
            end
    end
end

92 93 94 95 96 97 98 99 100 101 102
% 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')
    end
    if numel(zkl) ~= kl
        warning('Number of zkl values does not match the number specified in kl')
    end
end

103
if ftbverbose
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
    % 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)
    end

    % 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)
    end
123 124
end

125
% Calculate the sigma distributions at each grid node
126 127 128 129 130
switch lower(sigtype)
    case 'generalized'
        z = sigma_gen(nlev, dl, du, kl, ku, zkl, zku, ...
            Mobj.z(i), min_constant_depth);
    case 'uniform'
131 132 133 134
        z = sigma_geo(nlev, 1);
        % z = 0:-1/double(nlev-1):-1;
    case 'geometric'
        z = sigma_geo(nlev, sigpow);
135
    otherwise
136
        error('Don''t know sigtype %s (is it supported?)', sigtype)
137 138
end

139 140 141 142 143
% Create a siglay variable (i.e. midpoint in the sigma levels).
zlay = z(1:end-1) + (diff(z)/2);

Mobj.siglevz = repmat(Mobj.h, 1, nlev) .* repmat(z, Mobj.nVerts, 1);
Mobj.siglayz = repmat(Mobj.h, 1, nlev-1) .* repmat(zlay, Mobj.nVerts, 1);
144

145 146 147 148
% Add the sigma levels and layers to the Mobj.
Mobj.siglev = z;
Mobj.siglay = zlay;

149 150
if ftbverbose;
    fprintf(['end   : ' subname '\n'])
151
end