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