hybrid_coordinate

 Create a hybrid vertical coordinates file.

 Mobj = hybrid_coordinate(conf, Mobj);

 DESCRIPTION:
   Calculates and writes a hybird coordinates file for FVCOM.

 INPUT:
   conf - configuration struct with the following fields:
       sigma_file - file path into which to write the hybrid coordinates
       H0 - transition depth of the hybrid coordinates
       DU - upper water boundary thickness (metres)
       DL - lower water boundary thickness (metres)
       KU - number of layers in the DU water column
       KL - number of layers in the DL water column
       nlev - number of vertical levels (layers + 1)
   Mobj - Mesh object with the following fields:
       h - water depth at the nodes

 OUTPUT:
   Mobj - Mesh object with the following new fields:
       siglev - Sigma levels at the nodes
       siglevc - Sigma levels at the elements
       siglay - Sigma layers at the nodes
       siglayc - Sigma layers at the elements
       siglevz - Water depth levels at the nodes
       siglevzc - Water depth levels at the elements
       siglayz - Water depth layers at the nodes
       siglayzc - Water depth layers at the elements
       hc - Water depth at the elements

 EXAMPLE USAGE:
   conf.sigma_file = 'coord_hybrid.sig';
   conf.nlev = 41;
   conf.DU = 25;
   conf.DL = 25;
   conf.Hmin = 200;
   conf.KU = 5;
   conf.KL = 5;
   conf.ZKU = [.5 .5 .5 .5 .5];
   conf.ZKL = [.5 .5 .5 .5 .5];
   conf.H0 = 100;
   conf.nlev = 20;
   conf.DU = 25;
   conf.DL = 25;
   conf.KU = 5;
   conf.KL = 5;
   Mobj.h = random(100, 1) * 100;  % 100 random bathymetry points
   Mobj = hybrid_coordinate(conf, Mobj);

 Author(s):
   Ricard Torres (Plymouth Marine Laboratory)
   Pierre Cazenave (Plymouth Marine Laboratory)

 Revision history:
   2015-05-24 First version based on Riqui's initial implementation.
   2016-08-10 Updated the minimisation function to use the maximum of the
   difference between the two sets of vertical distributions rather than
   the median difference. Also tidy up the debug function.
   2017-01-26 Fix the transition depth optimisation and report the maximum
   difference between the two sigma level regions.

==========================================================================

0001 function Mobj = hybrid_coordinate(conf, Mobj)
0002 % Create a hybrid vertical coordinates file.
0003 %
0004 % Mobj = hybrid_coordinate(conf, Mobj);
0005 %
0006 % DESCRIPTION:
0007 %   Calculates and writes a hybird coordinates file for FVCOM.
0008 %
0009 % INPUT:
0010 %   conf - configuration struct with the following fields:
0011 %       sigma_file - file path into which to write the hybrid coordinates
0012 %       H0 - transition depth of the hybrid coordinates
0013 %       DU - upper water boundary thickness (metres)
0014 %       DL - lower water boundary thickness (metres)
0015 %       KU - number of layers in the DU water column
0016 %       KL - number of layers in the DL water column
0017 %       nlev - number of vertical levels (layers + 1)
0018 %   Mobj - Mesh object with the following fields:
0019 %       h - water depth at the nodes
0020 %
0021 % OUTPUT:
0022 %   Mobj - Mesh object with the following new fields:
0023 %       siglev - Sigma levels at the nodes
0024 %       siglevc - Sigma levels at the elements
0025 %       siglay - Sigma layers at the nodes
0026 %       siglayc - Sigma layers at the elements
0027 %       siglevz - Water depth levels at the nodes
0028 %       siglevzc - Water depth levels at the elements
0029 %       siglayz - Water depth layers at the nodes
0030 %       siglayzc - Water depth layers at the elements
0031 %       hc - Water depth at the elements
0032 %
0033 % EXAMPLE USAGE:
0034 %   conf.sigma_file = 'coord_hybrid.sig';
0035 %   conf.nlev = 41;
0036 %   conf.DU = 25;
0037 %   conf.DL = 25;
0038 %   conf.Hmin = 200;
0039 %   conf.KU = 5;
0040 %   conf.KL = 5;
0041 %   conf.ZKU = [.5 .5 .5 .5 .5];
0042 %   conf.ZKL = [.5 .5 .5 .5 .5];
0043 %   conf.H0 = 100;
0044 %   conf.nlev = 20;
0045 %   conf.DU = 25;
0046 %   conf.DL = 25;
0047 %   conf.KU = 5;
0048 %   conf.KL = 5;
0049 %   Mobj.h = random(100, 1) * 100;  % 100 random bathymetry points
0050 %   Mobj = hybrid_coordinate(conf, Mobj);
0051 %
0052 % Author(s):
0053 %   Ricard Torres (Plymouth Marine Laboratory)
0054 %   Pierre Cazenave (Plymouth Marine Laboratory)
0055 %
0056 % Revision history:
0057 %   2015-05-24 First version based on Riqui's initial implementation.
0058 %   2016-08-10 Updated the minimisation function to use the maximum of the
0059 %   difference between the two sets of vertical distributions rather than
0060 %   the median difference. Also tidy up the debug function.
0061 %   2017-01-26 Fix the transition depth optimisation and report the maximum
0062 %   difference between the two sigma level regions.
0063 %
0064 %==========================================================================
0065 
0066 [~, subname] = fileparts(mfilename('fullpath'));
0067 global ftbverbose
0068 if ftbverbose
0069     fprintf('\nbegin : %s\n', subname)
0070 end
0071 
0072 % Limits on the optimisation run.
0073 optimisation_settings = optimset('MaxFunEvals', 5000, ...
0074     'MaxIter', 5000, ...
0075     'TolFun', 10e-5, ...
0076     'TolX', 1e-7);
0077 
0078 % Extract the relevant information from the conf struct.
0079 nlev = conf.nlev;
0080 H0 = conf.H0;
0081 DU = conf.DU;
0082 DL = conf.DL;
0083 KU = conf.KU;
0084 KL = conf.KL;
0085 
0086 % Solve for Z0-Z2 to find Hmin parameter
0087 if ftbverbose
0088     fprintf('Optimising the hybrid coordinates... ')
0089 end
0090 ZKU = repmat(DU./KU, 1, KU);
0091 ZKL = repmat(DL./KL, 1, KL);
0092 fparams = @(H)hybrid_coordinate_hmin(H, nlev, DU, DL, KU, KL, ZKU, ZKL);
0093 [Hmin, minError] = fminsearch(fparams, H0, optimisation_settings);
0094 if ftbverbose
0095     fprintf('Hmin found %f with a max error in Vertical distribution of %f metres, \n',Hmin,minError)
0096     fprintf('Saving to %s... ', conf.sigma_file)
0097 end
0098 
0099 % Save to the given file name.
0100 fout = fopen(conf.sigma_file, 'wt');
0101 assert(fout >= 0, 'Error opening sigma file: %s', conf.sigma_file)
0102 fprintf(fout, 'NUMBER OF SIGMA LEVELS = %d\n', nlev);
0103 fprintf(fout, 'SIGMA COORDINATE TYPE = GENERALIZED\n');
0104 fprintf(fout, 'DU = %4.1f\n', DU);
0105 fprintf(fout, 'DL = %4.1f\n', DL);
0106 fprintf(fout, 'MIN CONSTANT DEPTH = %10.1f\n', round(Hmin));
0107 fprintf(fout, 'KU = %d\n', KU);
0108 fprintf(fout, 'KL = %d\n', KL);
0109 % Add the thicknesses with a loop.
0110 fprintf(fout, 'ZKU = ');
0111 for ii = 1:length(ZKU)
0112     fprintf(fout, '%4.1f', ZKU(ii));
0113 end
0114 fprintf(fout, '\n');
0115 fprintf(fout, 'ZKL = ');
0116 for ii = 1:length(ZKL)
0117     fprintf(fout, '%4.1f', ZKL(ii));
0118 end
0119 fprintf(fout,'\n');
0120 fclose(fout);
0121 
0122 if ftbverbose
0123     fprintf('done.\n')
0124     fprintf('Populating Mobj... ')
0125 end
0126 
0127 Mobj.siglev = zeros(Mobj.nVerts,nlev);
0128 Mobj.siglevc = zeros(Mobj.nElems,nlev);
0129 Mobj.siglayc = zeros(Mobj.nElems,nlev-1);
0130 
0131 % Create the arrays of the layer and level sigma coordinates.
0132 for xx = 1:length(Mobj.h)
0133     Mobj.siglev(xx,:) = sigma_gen(nlev,DL,DU,KL,KU,ZKL,ZKU,Mobj.h(xx),Hmin);
0134 end
0135 Mobj.siglay = Mobj.siglev(:,1:end-1) + (diff(Mobj.siglev,1,2)/2);
0136 % Turn off ftbverbose for this loop.
0137 old = ftbverbose;
0138 ftbverbose = 0;
0139 for zz = 1:nlev-1
0140     Mobj.siglevc(:, zz) = nodes2elems(Mobj.siglev(:, zz), Mobj);
0141     Mobj.siglayc(:, zz) = nodes2elems(Mobj.siglay(:, zz), Mobj);
0142 end
0143 % An extra level compared with layers.
0144 Mobj.siglevc(:, zz + 1) = nodes2elems(Mobj.siglev(:, zz + 1), Mobj);
0145 ftbverbose = old;
0146 clear old
0147 
0148 % Create a depth array for the element centres.
0149 if ~isfield(Mobj, 'hc')
0150     Mobj.hc = nodes2elems(Mobj.h, Mobj);
0151 end
0152 
0153 % Finally, make some [depth, sigma] arrays.
0154 Mobj.siglevz = repmat(Mobj.h, 1, nlev) .* Mobj.siglev;
0155 Mobj.siglayz = repmat(Mobj.h, 1, nlev-1) .* Mobj.siglay;
0156 Mobj.siglevzc = repmat(Mobj.hc, 1, nlev) .* Mobj.siglevc;
0157 Mobj.siglayzc = repmat(Mobj.hc, 1, nlev-1) .* Mobj.siglayc;
0158 
0159 if ftbverbose
0160     fprintf('done.\n')
0161     fprintf('end   : %s\n', subname)
0162 end
0163 
0164 return
0165 
0166 function ZZ = hybrid_coordinate_hmin(H, nlev, DU, DL, KU, KL, ZKU, ZKL)
0167 % Helper function to find the relevant minimum depth. I think.
0168 %
0169 %   ZZ = hybrid_coordinate_hmin(H, nlev, DU, DL, KU, KL, ZKU, ZKL)
0170 %
0171 % INPUT:
0172 %   H: transition depth of the hybrid coordinates?
0173 %   nlev: number of vertical levels (layers + 1)
0174 %   DU: upper water boundary thickness (metres)
0175 %   DL: lower water boundary thickness (metres)
0176 %   KU: layer number in the water column of DU
0177 %   KL: layer number in the water column of DL
0178 %
0179 % OUTPUT:
0180 %   ZZ: minimum water depth?
0181 %
0182 % Author(s):
0183 %   Ricard Torres (Plymouth Marine Laboratory)
0184 
0185 % if DU + DL > 1.25 * H;
0186 %     error('Depth %f too shallow for the chosen DU %f and DL %f values',H,DU,DL)
0187 % end
0188 
0189 Z0 = zeros(nlev, 1);
0190 Z2 = Z0;
0191 Z0(1, 1) = 0;
0192 DL2 = 0.001;
0193 DU2 = 0.001;
0194 KBM1 = nlev - 1;
0195 for nn = 1:nlev - 1
0196     X1 = DL2 + DU2;
0197     X1 = X1 * (KBM1 - nn) / KBM1;
0198     X1 = X1 - DL2;
0199     X1 = tanh(X1);
0200     X2 = tanh(DL2);
0201     X3 = X2 + tanh(DU2);
0202 
0203     Z0(nn + 1, 1)=((X1 + X2) / X3) - 1;
0204 end
0205 
0206 % s-coordinates
0207 X1 = (H - DU - DL);
0208 X2 = X1 ./ H;
0209 DR = X2 ./ (nlev - KU - KL - 1);
0210 
0211 Z2(1,1) = 0.0;
0212 
0213 for K = 2:KU + 1
0214     Z2(K, 1) = Z2(K - 1, 1) - (ZKU(K - 1) ./ H);
0215 end
0216 
0217 for K= KU + 2:nlev - KL
0218     Z2(K, 1) = Z2(K - 1, 1) - DR;
0219 end
0220 
0221 KK = 0;
0222 for K = nlev - KL + 1:nlev
0223     KK = KK + 1;
0224     Z2(K, 1) = Z2(K - 1, 1) - (ZKL(KK) ./ H);
0225 end
0226 ZZ = max(H*(Z0) - H*(Z2));
0227 % ZZ = max(abs(diff(H*(Z0)) - diff(H*(Z2))));
0228 
0229 return
0230 
0231 function debug_mode()
0232 % Test with made up data. This isn't actually used at all, but it's handy
0233 % to leave around for debugging things.
0234 
0235 conf.nlev = 25; % vertical levels (layers + 1)
0236 conf.H0 = 30; % threshold depth for the transition (metres)
0237 conf.DU = 3; % upper water boundary thickness
0238 conf.DL = 3; % lower water boundary thickness
0239 conf.KU = 3; % layer number in the water column of DU (maximum of 5 m thickness)
0240 conf.KL = 3; % layer number in the water column of DL (maximum of 5m thickness)
0241 
0242 
0243 Mobj = hybrid_coordinate(conf, Mobj);
0244 
0245 nlev = conf.nlev;
0246 H0 = conf.H0;
0247 DU = conf.DU;
0248 DL = conf.DL;
0249 KU = conf.KU;
0250 KL = conf.KL;
0251 ZKU = repmat(DU./KU, 1, KU);
0252 ZKL = repmat(DL./KL, 1, KL);
0253 
0254  Hmin=24;
0255 Hmax=Hmin + 200;
0256 y = 0:0.1:100;
0257 B = 70;
0258 H = Hmax .* exp(-((y./B-0.15).^2./0.5.^2));
0259 % H = [Hmin,H]; H=sort(H);
0260 nlev = conf.nlev;
0261 Z2=[];
0262 % Loop through all nodes to create sigma coordinates.
0263 for xx=1:length(H)
0264     Z2(xx, :) = sigma_gen(nlev, DL, DU, KL, KU, ZKL, ZKU, H(xx), Hmin);
0265 end
0266 
0267 clf
0268 plot(y,Z2 .* repmat(H', 1, nlev));hold on
0269 plot(y,ones(size(y)).*-Hmin)
0270 fprintf('Calculated minimum depth: %.2f\n', Hmin)
0271