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