0001 function [k,kc,Eyc,Q] = ADCP_DispCoef(beddepth,travdist,vertdepth,downstvel,startDist,endDist,extrp,extend_to_banks,banktype);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 nprof = length(beddepth);
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041 indx1 = find(isnan(beddepth));
0042 indx2 = find(~isnan(beddepth));
0043 beddepth(indx1) = interp1(travdist(indx2),beddepth(indx2),travdist(indx1));
0044
0045
0046
0047
0048
0049
0050
0051 for i = 1:nprof
0052 indx = find(~isnan(downstvel(i,:)));
0053 if isempty(indx)
0054 if i == 1 | i == nprof
0055 beep
0056 error('Edge NaN Detected in Downstream Velocity')
0057 end
0058 ustar(i) = nan;
0059 zo(i) = nan;
0060 else
0061 z = beddepth(i)-vertdepth(indx);
0062 zmean(i) = nanmean(z);
0063 u = downstvel(i,indx);
0064
0065
0066 [ustar1(i),zo(i)] = fitLogLawV2(u',z,beddepth(i));
0067 end
0068 end
0069
0070 ustar1 = real(ustar1);
0071 indx = find(isinf(ustar1));
0072 ustar1(indx) = nan;
0073
0074 indx1 = find(isnan(ustar1));
0075 indx2 = find(~isnan(ustar1));
0076 if isnan(ustar1(1))
0077 ustar1(1) = ustar1(indx2(1));
0078 end
0079 if isnan(ustar1(end))
0080 ustar1(end) = ustar1(indx2(end));
0081 end
0082 indx1 = find(isnan(ustar1));
0083 indx2 = find(~isnan(ustar1));
0084 ustar1(indx1) = interp1(travdist(indx2),ustar1(indx2),travdist(indx1));
0085
0086 zo = real(zo);
0087 indx = find(isinf(zo));
0088 zo(indx) = nan;
0089
0090 indx1 = find(isnan(zo));
0091 indx2 = find(~isnan(zo));
0092 if isnan(zo(1))
0093 zo(1) = zo(indx2(1));
0094 end
0095 if isnan(zo(end))
0096 zo(end) = zo(indx2(end));
0097 end
0098 indx1 = find(isnan(zo));
0099 indx2 = find(~isnan(zo));
0100 zo(indx1) = interp1(travdist(indx2),zo(indx2),travdist(indx1));
0101
0102 indx = find(zmean == 0);
0103 zmean(indx) = nan;
0104
0105 indx1 = find(isnan(zmean));
0106 indx2 = find(~isnan(zmean));
0107 if isnan(zmean(1))
0108 zmean(1) = zmean(indx2(1));
0109 end
0110 if isnan(zmean(end))
0111 zmean(end) = zmean(indx2(end));
0112 end
0113 indx1 = find(isnan(zmean));
0114 indx2 = find(~isnan(zmean));
0115 zmean(indx1) = interp1(travdist(indx2),zmean(indx2),travdist(indx1));
0116
0117 mzo = nanmedian(zo)
0118
0119
0120
0121
0122 test = downstvel;
0123 if extrp
0124 zgridspc = nanmean(diff(vertdepth));
0125 topDist = vertdepth(1)
0126 ntop = floor(topDist/zgridspc);
0127 top_nodes = linspace(0,topDist-zgridspc,ntop)';
0128 vertdepth = [top_nodes; vertdepth];
0129 downstvel = [nan*ones(length(beddepth),ntop) downstvel];
0130
0131
0132 for i = 1:nprof
0133 indx = find(~isnan(downstvel(i,:)));
0134
0135
0136 if ~isempty(indx)
0137
0138 top_nodes_flipped = beddepth(i) - top_nodes;
0139 u_top = ustar1(i)/0.41*log(top_nodes_flipped/zo(i));
0140
0141
0142 last_node = find(vertdepth < beddepth(i));
0143 bot_nodes = vertdepth(indx(end)+1:last_node(end));
0144
0145
0146
0147
0148
0149
0150
0151
0152 bot_nodes_flipped = beddepth(i) - bot_nodes;
0153 u_bot = ustar1(i)/0.41*log(bot_nodes_flipped/zo(i));
0154
0155
0156
0157
0158
0159 if ~isempty(indx)
0160 downstvel(i,indx(end)+1:last_node(end)) = u_bot;
0161 end
0162
0163
0164
0165 downstvel(i,1:ntop) = u_top;
0166
0167 end
0168
0169
0170
0171 clear indx last_node bot_nodes top_nodes_flipped bot_nodes_flipped u_top u_bot
0172 end
0173 end
0174
0175
0176
0177
0178
0179
0180
0181
0182 for i = 1:nprof
0183 lam(i) = VMT_LayerAveMean(vertdepth,downstvel(i,:)');
0184 end
0185
0186
0187
0188
0189 lam = real(lam)';
0190 indx = find(isinf(lam));
0191 lam(indx) = nan;
0192 indx1 = find(isnan(lam));
0193 indx2 = find(~isnan(lam));
0194 if isnan(lam(1))
0195 lam(1) = lam(indx2(1))/2;
0196 end
0197 if isnan(lam(end))
0198 lam(end) = lam(indx2(end))/2;
0199 end
0200 indx1 = find(isnan(lam));
0201 indx2 = find(~isnan(lam));
0202 lam(indx1) = interp1(travdist(indx2),lam(indx2),travdist(indx1));
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225 if 0
0226 shearvel = lam*0.41./log((beddepth/exp(1))./zoXS);
0227 else
0228 shearvel = lam*0.41./log((beddepth/exp(1))./mzo);
0229 end
0230
0231 indx = find(shearvel == 0);
0232 shearvel(indx) = nan;
0233
0234 indx1 = find(isnan(shearvel));
0235 indx2 = find(~isnan(shearvel));
0236 shearvel(indx1) = interp1(travdist(indx2),shearvel(indx2),travdist(indx1));
0237
0238
0239 travdist_orig = travdist;
0240
0241
0242 if extend_to_banks
0243 ygridspc = nanmean(diff(travdist));
0244 nstart = floor(startDist/ygridspc)
0245 nend = floor(endDist/ygridspc);
0246
0247
0248
0249 start_nodes = linspace(0,startDist-ygridspc,nstart);
0250 end_nodes = linspace(travdist(end)+ygridspc,travdist(end)+endDist,nend) + startDist;
0251
0252 switch banktype
0253 case 'triangular'
0254 startBed = interp1([start_nodes(1) startDist],[0 beddepth(1)],start_nodes);
0255 endBed = interp1([(travdist(end)+startDist) end_nodes(end)],[beddepth(end) 0],end_nodes);
0256 startLAM = lam(1).*sqrt(startBed./beddepth(1));
0257 endLAM = lam(end).*sqrt(endBed./beddepth(end));
0258 startSHV = interp1([start_nodes(1) startDist],[0 shearvel(1)],start_nodes);
0259 endSHV = interp1([(travdist(end)+startDist) end_nodes(end)],[shearvel(end) 0],end_nodes);
0260
0261
0262 case 'square'
0263 errordlg('Not implemented yet')
0264 end
0265
0266 travdist = [start_nodes'; travdist+startDist; end_nodes'];
0267 lam = [startLAM'; lam; endLAM'];
0268 beddepth = [startBed'; beddepth; endBed'];
0269 shearvel = [startSHV'; shearvel; endSHV'];
0270
0271 shearvel(1) = 0.1*shearvel(2);
0272 shearvel(end) = 0.1*shearvel(end-1);
0273
0274
0275
0276
0277
0278 end
0279
0280
0281 XSarea = trapz(travdist,beddepth)
0282
0283
0284
0285
0286 difftrav = diff(travdist);
0287 interim1 = 0.5*difftrav(1:end-1) + 0.5*difftrav(2:end);
0288 interim1 = [0.5*difftrav(1); interim1];
0289 bin_mp = cumsum(interim1);
0290
0291 interim2 = [0; bin_mp; travdist(end)];
0292 delW = diff(interim2);
0293
0294
0295
0296
0297 beddepth_mp = interp1(travdist,beddepth,bin_mp);
0298
0299
0300
0301
0302 cumlarea = cumtrapz(interim2,[beddepth(1); beddepth_mp; beddepth(end)]);
0303 bin_areas = diff(cumlarea);
0304
0305 indx1 = find(isnan(bin_areas));
0306 indx2 = find(~isnan(bin_areas));
0307 bin_areas(indx1) = interp1(travdist(indx2),bin_areas(indx2),travdist(indx1));
0308
0309
0310
0311
0312
0313 areaClosure = 100*(XSarea - sum(bin_areas))/XSarea
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360 d = bin_areas./delW;
0361
0362 indx1 = find(isnan(d));
0363 indx2 = find(~isnan(d));
0364 d(indx1) = interp1(travdist(indx2),d(indx2),travdist(indx1));
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376 Q = nansum(lam.*bin_areas)
0377 if Q < 0
0378 Q = -Q;
0379 lam = -lam;
0380 shearvel = -shearvel;
0381 end
0382
0383
0384 Ua = Q./XSarea
0385
0386
0387 uprime = lam - Ua;
0388
0389
0390
0391 B = travdist(end) - travdist(1)
0392 if 1
0393 Ey = 0.21*shearvel.*d;
0394
0395 indx1 = find(isnan(Ey));
0396 indx2 = find(~isnan(Ey));
0397 Ey(indx1) = interp1(travdist(indx2),Ey(indx2),travdist(indx1));
0398 else
0399
0400 Ustar = trapz(travdist,shearvel)/B
0401 H = XSarea/B
0402
0403 theta = 0.145 + (1/3520)*(Ua/Ustar)*(B/H)^1.38
0404
0405 Ey = theta*Ustar*d;
0406
0407
0408 indx1 = find(isnan(Ey));
0409 indx2 = find(~isnan(Ey));
0410 Ey(indx1) = interp1(travdist(indx2),Ey(indx2),travdist(indx1));
0411 end
0412
0413 Eyc = trapz(travdist,Ey)/B;
0414
0415
0416
0417
0418
0419 INT1 = cumtrapz(travdist,uprime.*d);
0420
0421
0422
0423 INT2 = cumtrapz(travdist,INT1./(Ey.*d));
0424
0425
0426
0427 INT3 = cumtrapz(travdist,INT2.*uprime.*d);
0428
0429
0430 if 0
0431 figure(10); clf;
0432 title('Variable Ey')
0433 subplot(3,1,1)
0434 plot(travdist,INT1,'b-')
0435 ylabel('INT1')
0436 subplot(3,1,2); plot(travdist,INT2,'b-')
0437 ylabel('INT2')
0438 subplot(3,1,3); plot(travdist,INT3,'b-')
0439 ylabel('INT3')
0440 end
0441
0442
0443
0444 k = -INT3(end)/XSarea;
0445
0446 disp(' ')
0447 disp(['The longitudinal dispersion coefficient(K)is ' num2str(k) ' m^2/s'])
0448 disp(' with a VARIABLE transverse mixing coefficient(Ey) ')
0449
0450
0451
0452
0453
0454
0455 INT2c = cumtrapz(travdist,INT1./(Eyc.*d));
0456
0457
0458
0459 INT3c = cumtrapz(travdist,INT2c.*uprime.*d);
0460
0461
0462 if 0
0463 figure(11); clf; title('Constant Ey')
0464 subplot(2,1,1); plot(travdist,INT2c,'b-')
0465 ylabel('INT2')
0466 subplot(2,1,2); plot(travdist,INT3c,'b-')
0467 ylabel('INT3')
0468 end
0469
0470
0471
0472 kc = -INT3c(end)/XSarea;
0473
0474 disp(' ')
0475 disp(['The longitudinal dispersion coefficient(K)is ' num2str(kc) ' m^2/s'])
0476 disp(' with a CONSTANT transverse mixing coefficient(Ey) ')
0477
0478
0479 figure(8); clf
0480 subplot(3,1,1); plot(travdist,beddepth,'k-',travdist,d,'r-'); ylabel('Bed Depth (m)')
0481
0482 subplot(3,1,2); plot(travdist,lam,'k-'); ylabel('Depth-Averaged Vel. (m/s)')
0483 subplot(3,1,3); plot(travdist,shearvel,'k-'); ylabel('u* (m/s)')
0484 xlabel('Transverse Distance (m)')
0485 figure(9); clf
0486 subplot(3,1,1); plot(travdist,bin_areas,'k-'); ylabel('Bin Areas (m^2)')
0487 subplot(3,1,2); plot(travdist,uprime,'k-'); ylabel('uprime (m/s)')
0488 subplot(3,1,3); plot(travdist,Ey,'k-'); ylabel('Ey (m^2/s)')
0489 xlabel('Transverse Distance (m)')
0490
0491
0492
0493
0494
0495
0496
0497