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.

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

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 %
0062 %==========================================================================
0063 
0064 [~, subname] = fileparts(mfilename('fullpath'));
0065 global ftbverbose
0066 if ftbverbose
0067     fprintf('\nbegin : %s\n', subname)
0068 end
0069 
0070 % Limits on the optimisation run.
0071 optimisation_settings = optimset('MaxFunEvals', 5000, ...
0072     'MaxIter', 5000, ...
0073     'TolFun', 10e-5, ...
0074     'TolX', 1e-7);
0075 
0076 % Extract the relevant information from the conf struct.
0077 nlev = conf.nlev;
0078 H0 = conf.H0;
0079 DU = conf.DU;
0080 DL = conf.DL;
0081 KU = conf.KU;
0082 KL = conf.KL;
0083 
0084 % Solve for Z0-Z2 to find Hmin parameter
0085 if ftbverbose
0086     fprintf('Optimising the hybrid coordinates... ')
0087 end
0088 ZKU = repmat(DU./KU, 1, KU);
0089 ZKL = repmat(DL./KL, 1, KL);
0090 fparams = @(H)hybrid_coordinate_hmin(H, nlev, DU, DL, KU, KL, ZKU, ZKL);
0091 [Hmin, ~] = fminsearch(fparams, H0, optimisation_settings);
0092 
0093 if ftbverbose
0094     fprintf('done.\n')
0095     fprintf('Saving to %s... ', conf.sigma_file)
0096 end
0097 
0098 % Save to the given file name.
0099 fout = fopen(conf.sigma_file, 'wt');
0100 assert(fout >= 0, 'Error opening sigma file: %s', conf.sigma_file)
0101 fprintf(fout, 'NUMBER OF SIGMA LEVELS = %d\n', nlev);
0102 fprintf(fout, 'SIGMA COORDINATE TYPE = GENERALIZED\n');
0103 fprintf(fout, 'DU = %4.1f\n', DU);
0104 fprintf(fout, 'DL = %4.1f\n', DL);
0105 fprintf(fout, 'MIN CONSTANT DEPTH = %10.1f\n', round(Hmin));
0106 fprintf(fout, 'KU = %d\n', KU);
0107 fprintf(fout, 'KL = %d\n', KL);
0108 % Add the thicknesses with a loop.
0109 fprintf(fout, 'ZKU = ');
0110 for ii = 1:length(ZKU)
0111     fprintf(fout, '%4.1f', ZKU(ii));
0112 end
0113 fprintf(fout, '\n');
0114 fprintf(fout, 'ZKL = ');
0115 for ii = 1:length(ZKU)
0116     fprintf(fout, '%4.1f', ZKL(ii));
0117 end
0118 fprintf(fout,'\n');
0119 fclose(fout);
0120 
0121 if ftbverbose
0122     fprintf('done.\n')
0123     fprintf('Populating Mobj... ')
0124 end
0125 
0126 % Create the arrays of the layer and level sigma coordinates.
0127 for xx = 1:length(Mobj.h)
0128     Mobj.siglev(xx,:) = sigma_gen(nlev,DL,DU,KL,KU,ZKL,ZKU,Mobj.h(xx),Hmin);
0129 end
0130 Mobj.siglay = Mobj.siglev(:,1:end-1) + (diff(Mobj.siglev,1,2)/2);
0131 % Turn off ftbverbose for this loop.
0132 old = ftbverbose;
0133 ftbverbose = 0;
0134 for zz = 1:nlev-1
0135     Mobj.siglevc(:, zz) = nodes2elems(Mobj.siglev(:, zz), Mobj);
0136     Mobj.siglayc(:, zz) = nodes2elems(Mobj.siglay(:, zz), Mobj);
0137 end
0138 % An extra level compared with layers.
0139 Mobj.siglevc(:, zz + 1) = nodes2elems(Mobj.siglev(:, zz + 1), Mobj);
0140 ftbverbose = old;
0141 clear old
0142 
0143 % Create a depth array for the element centres.
0144 if ~isfield(Mobj, 'hc')
0145     Mobj.hc = nodes2elems(Mobj.h, Mobj);
0146 end
0147 
0148 % Finally, make some [depth, sigma] arrays.
0149 Mobj.siglevz = repmat(Mobj.h, 1, nlev) .* Mobj.siglev;
0150 Mobj.siglayz = repmat(Mobj.h, 1, nlev-1) .* Mobj.siglay;
0151 Mobj.siglevzc = repmat(Mobj.hc, 1, nlev) .* Mobj.siglevc;
0152 Mobj.siglayzc = repmat(Mobj.hc, 1, nlev-1) .* Mobj.siglayc;
0153 
0154 if ftbverbose
0155     fprintf('done.\n')
0156     fprintf('end   : %s\n', subname)
0157 end
0158 
0159 return
0160 
0161 function ZZ = hybrid_coordinate_hmin(H, nlev, DU, DL, KU, KL, ZKU, ZKL)
0162 % Helper function to find the relevant minimum depth. I think.
0163 %
0164 %   ZZ = hybrid_coordinate_hmin(H, nlev, DU, DL, KU, KL, ZKU, ZKL)
0165 %
0166 % INPUT:
0167 %   H: transition depth of the hybrid coordinates?
0168 %   nlev: number of vertical levels (layers + 1)
0169 %   DU: upper water boundary thickness (metres)
0170 %   DL: lower water boundary thickness (metres)
0171 %   KU: layer number in the water column of DU
0172 %   KL: layer number in the water column of DL
0173 %
0174 % OUTPUT:
0175 %   ZZ: minimum water depth?
0176 %
0177 % Author(s):
0178 %   Ricard Torres (Plymouth Marine Laboratory)
0179 
0180 if DU + DL > 1.25 * H;
0181     disp('Depth too shallow for the chosen DU and DL values')
0182     return
0183 end
0184 
0185 Z0 = zeros(nlev, 1);
0186 Z2 = Z0;
0187 Z0(1, 1) = 0;
0188 DL2 = 0.001;
0189 DU2 = 0.001;
0190 KBM1 = nlev - 1;
0191 for nn = 1:nlev - 1
0192     X1 = DL2 + DU2;
0193     X1 = X1 * (KBM1 - nn) / KBM1;
0194     X1 = X1 - DL2;
0195     X1 = tanh(X1);
0196     X2 = tanh(DL2);
0197     X3 = X2 + tanh(DU2);
0198 
0199     Z0(nn + 1, 1)=((X1 + X2) / X3) - 1;
0200 end
0201 
0202 % s-coordinates
0203 X1 = (H - DU - DL);
0204 X2 = X1 ./ H;
0205 DR = X2 ./ (nlev - KU - KL - 1);
0206 
0207 Z2(1,1) = 0.0;
0208 
0209 for K = 2:KU + 1
0210     Z2(K, 1) = Z2(K - 1, 1) - (ZKU(K - 1) ./ H);
0211 end
0212 
0213 for K= KU + 2:nlev - KL
0214     Z2(K, 1) = Z2(K - 1, 1) - DR;
0215 end
0216 
0217 KK = 0;
0218 for K = nlev - KL + 1:nlev
0219     KK = KK + 1;
0220     Z2(K, 1) = Z2(K - 1, 1) - (ZKL(KK) ./ H);
0221 end
0222 ZZ = max(diff(Z0) - diff(Z2));
0223 
0224 return
0225 
0226 function debug_mode()
0227 % Test with made up data. This isn't actually used at all, but it's handy
0228 % to leave around for debugging things.
0229 
0230 Hmin=50;
0231 Hmax=Hmin + 200;
0232 y = 0:0.1:100;
0233 B = 100;
0234 H = Hmax .* exp(-((y./B-0.15).^2./0.5.^2));
0235 nlev = conf.nlev;
0236 
0237 % Loop through all nodes to create sigma coordinates.
0238 for xx=1:length(H)
0239     Z2(xx, :) = sigma_gen(nlev, DL, DU, KL, KU, ZKL, ZKU, H(xx), Hmin);
0240 end
0241 
0242 plot(Z2 .* repmat(H', 1, nlev))
0243 fprintf('Calculated minimum depth: %.2f\n', Hmin)
0244