read_sigma

 Read an FVCOM sigma layers file and output z values into Mobj.
 
 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:
   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.
 
 EXAMPLE USAGE:
   Mobj = read_sigma(Mobj, 'sigma.dat')
 
 Author(s):
   Pierre Cazenave (Plymouth Marine Laboratory)
 
 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
   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.

0001 function Mobj = read_sigma(Mobj, sigmafile)
0002 % Read an FVCOM sigma layers file and output z values into Mobj.
0003 %
0004 % Mobj = read_sigma(Mobj, sigmafile)
0005 %
0006 % DESCRIPTION:
0007 %   Read a sigma file and calculate the sigma layer depths
0008 %
0009 % INPUT:
0010 %   Mobj:       Mesh object which must contain Mobj.h (depths).
0011 %   sigmafile : Full path to an FVCOM sigma.dat file.
0012 %
0013 % OUTPUT:
0014 %   Mobj:       Mesh object with four new fields:
0015 %               - siglayz and siglevz: contain depths of the sigma layers
0016 %               and levels at each grid node.
0017 %               - siglay and siglev: the sigma layer and levels in the
0018 %               range 0 to -1.
0019 %
0020 % EXAMPLE USAGE:
0021 %   Mobj = read_sigma(Mobj, 'sigma.dat')
0022 %
0023 % Author(s):
0024 %   Pierre Cazenave (Plymouth Marine Laboratory)
0025 %
0026 % Revision history
0027 %   2013-01-08 Based on the code in show_sigma.m but instead of calculating
0028 %   sigma layers along a user-defined line, the depths are calculated for
0029 %   each node in the unstructured grid.
0030 %   2013-01-10 Added two new outputs to the returned Mobj (siglay and
0031 %   siglev). They're useful in write_FVCOM_tsobc.m.
0032 
0033 subname = 'read_sigma';
0034 
0035 global ftbverbose;
0036 if ftbverbose
0037     fprintf('\n')
0038     fprintf(['begin : ' subname '\n'])
0039 end
0040 
0041 fid = fopen(sigmafile,'r');
0042 if(fid  < 0)
0043     error(['file: ' sigmafile ' does not exist']);
0044 end
0045 
0046 while ~feof(fid)
0047     line = fgetl(fid);
0048     if isempty(line) || strncmp(line, '!', 1) || ~ischar(line)
0049         continue
0050     end
0051     key = lower(line(1:3));
0052     C = strtrim(regexpi(line, '=', 'split'));
0053     switch key
0054         case 'num'
0055             nlev = str2double(C{2});
0056         case 'sig'
0057             sigtype = C{2};
0058         case 'du '
0059             du = str2double(C{2});
0060         case 'dl '
0061             dl = str2double(C{2});
0062         case 'min'
0063             min_constant_depth = str2double(C{2});
0064         case 'ku '
0065             ku = str2double(C{2});
0066         case 'kl '
0067             kl = str2double(C{2});
0068         case 'zku'
0069             s = str2double(regexp(C{2}, ' ', 'split'));
0070             zku = zeros(ku, 1);
0071             for i = 1:ku
0072                 zku(i) = s(i);
0073             end
0074         case 'zkl'
0075             s = str2double(regexp(C{2}, ' ', 'split'));
0076             zkl = zeros(kl, 1);
0077             for i = 1:kl
0078                 zkl(i) = s(i);
0079             end
0080     end
0081 end
0082 
0083 if ftbverbose
0084     fprintf('nlev %d\n',nlev)
0085     fprintf('sigtype %s\n',sigtype)
0086     fprintf('du %d\n',du)
0087     fprintf('dl %d\n',dl)
0088     fprintf('min_constant_depth %f\n',min_constant_depth)
0089     fprintf('ku %d\n',ku)
0090     fprintf('kl %d\n',kl)
0091     fprintf('zku %d\n',zku)
0092     fprintf('zkl %d\n',zkl)
0093 end
0094 
0095 % calculate the sigma distributions at each grid node
0096 switch lower(sigtype)
0097     case 'generalized'
0098         z = sigma_gen(nlev, dl, du, kl, ku, zkl, zku, ...
0099             Mobj.z(i), min_constant_depth);
0100     case 'uniform'
0101         z = 0:-1/double(nlev-1):-1;
0102     otherwise
0103         error('Can''t do that sigtype')
0104 end
0105 
0106 % Create a siglay variable (i.e. midpoint in the sigma levels).
0107 zlay = z(1:end-1) + (diff(z)/2);
0108 
0109 Mobj.siglevz = repmat(Mobj.h, 1, nlev) .* repmat(z, Mobj.nVerts, 1);
0110 Mobj.siglayz = repmat(Mobj.h, 1, nlev-1) .* repmat(zlay, Mobj.nVerts, 1);
0111 
0112 % Add the sigma levels and layers to the Mobj.
0113 Mobj.siglev = z;
0114 Mobj.siglay = zlay;
0115 
0116 if ftbverbose;
0117     fprintf(['end   : ' subname '\n'])
0118 end