# ADCP_DispCoef

## PURPOSE

This program computes the longitudinal dispersion coefficient from ADCP

## SYNOPSIS

function [k,kc,Eyc,Q] = ADCP_DispCoef(beddepth,travdist,vertdepth,downstvel,startDist,endDist,extrp,extend_to_banks,banktype);

## DESCRIPTION

    This program computes the longitudinal dispersion coefficient from ADCP
    transects.

## CROSS-REFERENCE INFORMATION

This function calls:

- fitLogLawV2 This function fits a log law of the form u/u* = 1/kappa*ln(z/zo) to the given data
- nansum NANSUM Sum, ignoring NaNs.

This function is called by:

## SOURCE CODE

    0001 function [k,kc,Eyc,Q] = ADCP_DispCoef(beddepth,travdist,vertdepth,downstvel,startDist,endDist,extrp,extend_to_banks,banktype);
    0002
    0003 %This program computes the longitudinal dispersion coefficient from ADCP
    0004 %transects.
    0005
    0006 %P.R. Jackson & N.V. Reynolds, USGS, 11/16/10
    0007
    0008 %Inputs:
    0009 %
    0010 %     beddepth  = Depth (in meters) from the water surface to the bed (n element vector)
    0011 %     travdist  = Transverse distance (in m) across the river starting with (n element vector)
    0012 %     vertdepth = Depths (in m) from the water surface to velocity bins in the water column (m element vector)
    0013 %     downstvel = streamwise velocity distribution (+ for DS) in m/s (n x m array)
    0014 %     startDist = distance (in m) from the first profile to the starting bank
    0015 %     endDist   = distance (in m) from the last profile to the ending bank
    0016 %     extrp     = (Binary) Set to 1 to extrapolate profiles to the bed and surface (log law); otherwise set to 0
    0017 %     extend_to_banks  = (Binary) Set to 1 to extrapolate data to the banks; otherwise set to 0
    0018 %     banktype  = (string) 'Triangular' or 'Square' (used in bank extrapolation)
    0019
    0020 %Outputs:
    0021 %     k  = Longitudinal dispersion coefficient (in m^2/s) (variable Ey)
    0022 %     kc  = Longitudinal dispersion coefficient (in m^2/s) (Constant Ey)
    0023 %     Ey = Transverse mixing coefficient (constant)
    0024 %     Q  = Approximate discharge in m^3/s
    0025
    0026
    0027
    0028 %% Basic parameters
    0029 nprof = length(beddepth);
    0030
    0031 %% Correct missing bed depths
    0032 %figure(10); clf; plot(travdist,beddepth,'k-')
    0033 % indx1 = find(isnan(beddepth));
    0034 % indx2 = find(~isnan(beddepth));
    0035 % if isnan(beddepth(1))
    0036 %     beddepth(1) = beddepth(indx2(1))/2; %Fill end nans (not filled with interpolation)
    0037 % end
    0038 % if isnan(beddepth(end))
    0039 %     beddepth(end) = beddepth(indx2(end))/2;  %Fill end nans (not filled with interpolation)
    0040 % end
    0041 indx1 = find(isnan(beddepth));
    0042 indx2 = find(~isnan(beddepth));
    0043 beddepth(indx1) = interp1(travdist(indx2),beddepth(indx2),travdist(indx1));  %fills the nans
    0044
    0045 %figure(10); hold on; plot(travdist(indx1),beddepth(indx1),'r-')
    0046
    0047 %% Fit each profile with a log law
    0048
    0049 %[ustarXS,zoXS] = ADCP_ComputeXS_ShearVelocity(downstvel,vertdepth,beddepth)
    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);  % Height above the bed
    0062         zmean(i) = nanmean(z);  %Mean height above the bed
    0063         u = downstvel(i,indx);
    0064         %figure(10); clf; plot(downstvel(i,indx),beddepth(i)-vertdepth(indx))
    0065         %[ustar(i),zo(i),ks,cod,upred,zpred,delta] = fitLogLaw(u',z,beddepth(i));
    0066         [ustar1(i),zo(i)] = fitLogLawV2(u',z,beddepth(i));
    0067     end
    0068 end
    0069 %clean up problem values
    0070 ustar1 = real(ustar1);
    0071 indx = find(isinf(ustar1)); %turn if values to NAN
    0072 ustar1(indx) = nan;
    0073 %edge values
    0074 indx1 = find(isnan(ustar1));
    0075 indx2 = find(~isnan(ustar1));
    0076 if isnan(ustar1(1))
    0077     ustar1(1) = ustar1(indx2(1)); %Fill end nans (not filled with interpolation)
    0078 end
    0079 if isnan(ustar1(end))
    0080     ustar1(end) = ustar1(indx2(end));  %Fill end nans (not filled with interpolation)
    0081 end
    0082 indx1 = find(isnan(ustar1));
    0083 indx2 = find(~isnan(ustar1));
    0084 ustar1(indx1) = interp1(travdist(indx2),ustar1(indx2),travdist(indx1));  %fills the nans
    0085
    0086 zo = real(zo);
    0087 indx = find(isinf(zo)); %turn if values to NAN
    0088 zo(indx) = nan;
    0089 %edge values
    0090 indx1 = find(isnan(zo));
    0091 indx2 = find(~isnan(zo));
    0092 if isnan(zo(1))
    0093     zo(1) = zo(indx2(1)); %Fill end nans (not filled with interpolation)
    0094 end
    0095 if isnan(zo(end))
    0096     zo(end) = zo(indx2(end));  %Fill end nans (not filled with interpolation)
    0097 end
    0098 indx1 = find(isnan(zo));
    0099 indx2 = find(~isnan(zo));
    0100 zo(indx1) = interp1(travdist(indx2),zo(indx2),travdist(indx1));  %fills the nans
    0101
    0102 indx = find(zmean == 0); %turn if values to NAN
    0103 zmean(indx) = nan;
    0104 %edge values
    0105 indx1 = find(isnan(zmean));
    0106 indx2 = find(~isnan(zmean));
    0107 if isnan(zmean(1))
    0108     zmean(1) = zmean(indx2(1)); %Fill end nans (not filled with interpolation)
    0109 end
    0110 if isnan(zmean(end))
    0111     zmean(end) = zmean(indx2(end));  %Fill end nans (not filled with interpolation)
    0112 end
    0113 indx1 = find(isnan(zmean));
    0114 indx2 = find(~isnan(zmean));
    0115 zmean(indx1) = interp1(travdist(indx2),zmean(indx2),travdist(indx1));  %fills the nans
    0116
    0117 mzo = nanmedian(zo)  %Median z0
    0118 %figure(10); clf; plot(travdist,ustar','k',travdist,shearvel,'r',travdist,ustar2','b')
    0119
    0120 %% Extrapolate top and bottom portions if required
    0121 %figure(20); clf; contour(downstvel,'Fill','on','Linestyle','none');
    0122 test = downstvel;
    0123 if extrp  %extrapolates to the bed and surface
    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];  %preallocated space for new top velocities
    0130     %figure(21); clf; contour(downstvel,'Fill','on','Linestyle','none');
    0131
    0132     for i = 1:nprof
    0133         indx = find(~isnan(downstvel(i,:)));
    0134
    0135         %figure(10); clf; plot(downstvel(i,:),vertdepth,'ko-'); set(gca,'YDir','reverse');
    0136         if ~isempty(indx) %all nans get no computations or extensions
    0137
    0138             top_nodes_flipped = beddepth(i) - top_nodes;
    0139             u_top = ustar1(i)/0.41*log(top_nodes_flipped/zo(i));
    0140             %botDist = beddepth(i) - vertdepth(indx(end));
    0141
    0142             last_node = find(vertdepth < beddepth(i));
    0143             bot_nodes = vertdepth(indx(end)+1:last_node(end));
    0144
    0145             %nbot   = floor(botDist/zgridspc);
    0146             %bot_nodes = linspace(vertdepth(indx(end))+zgridspc,vertdepth(indx(end))+botDist,nbot);
    0147
    0148             %Flip the vertical coordinate and use log fit from above to extend
    0149             %to boundaries
    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             %u_bot(end) = 0.0;  %Log law is singular at z = 0, so reset to zero for no slip
    0156
    0157
    0158             % Append to data
    0159             if ~isempty(indx)
    0160                 downstvel(i,indx(end)+1:last_node(end)) = u_bot;
    0161             end
    0162             %downstvel(i,last_node(end)+1) = 0.0;  %Set to zero at bed
    0163             %vertdepth(last_node(end)+1) = beddepth(i); %reset first invalid
    0164             %bottom bin to bed depth
    0165             downstvel(i,1:ntop) = u_top;
    0166
    0167         end
    0168
    0169         %figure(10); hold on; plot(u_top,top_nodes,'r.',u_bot,bot_nodes,'b.'); set(gca,'YDir','reverse'); pause
    0170
    0171         clear indx last_node bot_nodes top_nodes_flipped bot_nodes_flipped u_top u_bot
    0172     end
    0173 end
    0174 %test2 = downstvel(:,ntop+1:end) - test;
    0175 %figure(22); clf; contour(test2,'Fill','on','Linestyle','none');
    0176
    0177 %% Compute Depth-averaged velocities
    0178
    0179 %This function computes the layer averaged mean of the downstream velocity over the measured depth range.
    0180 %Assumes the data outside the depth range have been set to NaN.
    0181 % P.R. Jackson, USGS 1-7-09
    0182 for i = 1:nprof
    0183     lam(i) = VMT_LayerAveMean(vertdepth,downstvel(i,:)');
    0184 end
    0185
    0186 %figure(10); clf; plot(travdist,lam,'r-') % Plot to check
    0187
    0188 %lam = VMT_LayerAveMean(repmat(vertdepth,1,size(downstvel',2)),downstvel')';
    0189 lam = real(lam)';  %remove any imaginary parts
    0190 indx = find(isinf(lam)); %turn if values to NAN
    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; %Fill end nans (not filled with interpolation)
    0196 end
    0197 if isnan(lam(end))
    0198     lam(end) = lam(indx2(end))/2;  %Fill end nans (not filled with interpolation)
    0199 end
    0200 indx1 = find(isnan(lam));
    0201 indx2 = find(~isnan(lam));
    0202 lam(indx1) = interp1(travdist(indx2),lam(indx2),travdist(indx1));  %fills the nans
    0203
    0204 %indx1 = find(isnan(lam))
    0205 %figure(10); clf; plot(travdist,lam,'r-') % Plot to check
    0206
    0207 % %Davide's method for LAM
    0208 %
    0209 % for i=1:sztravdist+1
    0210 %     x(i)=nansum(downstvel(i,:));
    0211 %
    0212 % end
    0213 %
    0214 %
    0215 % lamd=x'./beddepth;
    0216 % lamd(1)=lamd(2)/2;
    0217 % lamd(sztravdist+1)=lamd(sztravdist)/2;
    0218 % %lamd=lamd';
    0219 % figure(10); hold on; plot(travdist,lamd./3.281,'b-')
    0220 % whos
    0221 % return
    0222
    0223 %% Recompute the shear velocity using the median z0 and LAM (lower noise)
    0224 %shearvel = lam*0.41./log(zmean'./mzo);  %kappa = 0.41
    0225 if 0  %use the z0 from the normalized fit
    0226     shearvel = lam*0.41./log((beddepth/exp(1))./zoXS);  %kappa = 0.41; %Use zmean = h/e following Sime etal 2007 as mean velocity occurs at h/e
    0227 else %uses the median of the individual profile fits
    0228     shearvel = lam*0.41./log((beddepth/exp(1))./mzo);  %kappa = 0.41; %Use zmean = h/e following Sime etal 2007 as mean velocity occurs at h/e
    0229 end
    0230 % replace zeros with nans
    0231 indx = find(shearvel == 0);
    0232 shearvel(indx) = nan;
    0233 % Fill gaps
    0234 indx1 = find(isnan(shearvel));
    0235 indx2 = find(~isnan(shearvel));
    0236 shearvel(indx1) = interp1(travdist(indx2),shearvel(indx2),travdist(indx1));  %fills the nans
    0237
    0238 %figure(10); clf; plot(travdist,shearvel,'r')
    0239 travdist_orig = travdist;
    0240
    0241 %% Extrapolate to the banks (following WRII)
    0242 if extend_to_banks
    0243     ygridspc = nanmean(diff(travdist));
    0244     nstart = floor(startDist/ygridspc)
    0245     nend   = floor(endDist/ygridspc);
    0246 %Note: This fails for zero start or end distances
    0247     %indx = find(~isnan(lam));
    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);  %Linear interpolation of shear velocities
    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);  %Reset start and end u* to nonzero to keep INT2 from blowing up (due to zero Ey at edges)
    0272     shearvel(end) = 0.1*shearvel(end-1);
    0273
    0274
    0275     %figure(10); hold on; plot(travdist,lam,'b-') % Plot to check
    0276     %figure(10); clf; plot(travdist,beddepth,'b-') % Plot to check
    0277
    0278 end
    0279
    0280 %% Compute the mean Cross-sectional area
    0281 XSarea = trapz(travdist,beddepth)
    0282
    0283
    0284 %% Compute bin areas
    0285 %compute the midpoints
    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 %figure(10); clf; plot(travdist,beddepth,'k.-');
    0295
    0296 %Interpolate the bed depths at the midpoints
    0297 beddepth_mp = interp1(travdist,beddepth,bin_mp);
    0298
    0299 %figure(10); clf; plot(travdist,beddepth,'k.-',bin_mp,beddepth_mp,'r.'); %plot to check
    0300
    0301 %Compute the areas
    0302 cumlarea = cumtrapz(interim2,[beddepth(1); beddepth_mp; beddepth(end)]);
    0303 bin_areas = diff(cumlarea);
    0304 % Fill gaps
    0305 indx1 = find(isnan(bin_areas));
    0306 indx2 = find(~isnan(bin_areas));
    0307 bin_areas(indx1) = interp1(travdist(indx2),bin_areas(indx2),travdist(indx1));  %fills the nans
    0308
    0309
    0310 %figure(10); clf; plot(travdist,bin_areas,'k.-'); %plot to check
    0311
    0312 %Check the computation
    0313 areaClosure = 100*(XSarea - sum(bin_areas))/XSarea
    0314
    0315
    0316     %Old area comp
    0317     % sztravdist=size(travdist,1)-1;
    0318     % szvertdepth=size(vertdepth,1);
    0319     %
    0320     % %midpoint to the left of y-component
    0321     % midptly=(1:sztravdist+1);
    0322     % midptly(1)=0;
    0323     % midptly(2)=travdist(2)-(travdist(3)-travdist(2))/2;
    0324     % midptly(3:sztravdist-0)=(travdist(3:sztravdist)+travdist(2:sztravdist-1))/2;
    0325     % midptly(sztravdist+1)=travdist(sztravdist)-(travdist(sztravdist-1)-travdist(sztravdist))/2;
    0326     %
    0327     % %midpoint to the left of bed depth
    0328     % midptld=(1:sztravdist+1);
    0329     % midptld(1)=0;
    0330     % midptld(2)=beddepth(2)-(beddepth(3)-beddepth(2))/2;
    0331     % midptld(3:sztravdist-0)=(beddepth(3:sztravdist)+beddepth(2:sztravdist-1))/2;
    0332     % midptld(sztravdist+1)=beddepth(sztravdist)-(beddepth(sztravdist-1)-beddepth(sztravdist))/2;
    0333     %
    0334     % %midpoint to the right of y-component
    0335     % midptry=(1:sztravdist+1);
    0336     % midptry(1)=travdist(2)-(travdist(3)-travdist(2))/2;
    0337     % midptry(2:sztravdist-1)=(travdist(3:sztravdist)+travdist(2:sztravdist-1))/2;
    0338     % midptry(sztravdist)=travdist(sztravdist)-(travdist(sztravdist-1)-travdist(sztravdist))/2;
    0339     % midptry(sztravdist+1)=0;
    0340     %
    0341     % %midpoint to the right of bed depth
    0342     % midptrd=(1:sztravdist+1);
    0343     % midptrd(1)=beddepth(2)-(beddepth(3)-beddepth(2))/2;
    0344     % midptrd(2:sztravdist-1)=(beddepth(3:sztravdist)+beddepth(2:sztravdist-1))/2;
    0345     % midptrd(sztravdist)=beddepth(sztravdist)-(beddepth(sztravdist-1)-beddepth(sztravdist))/2;
    0346     % midptrd(sztravdist+1)=0;
    0347     %
    0348     %
    0349     % %area per bin
    0350     % ameth2=(1:sztravdist+1);
    0351     % ameth2(1)=(beddepth(1)+midptrd(1))*(midptry(1)-travdist(1))/2;
    0352     % for i=2:sztravdist
    0353     % ameth2(i)=(beddepth(i)+midptld(i)')*(travdist(i)-midptly(i)')/2+(beddepth(i)+midptrd(i)').*(midptry(i)'-travdist(i))./2;
    0354     % end
    0355     % ameth2(sztravdist+1)=((beddepth(sztravdist+1)+midptld(sztravdist+1))*(travdist(sztravdist+1)-midptly(sztravdist+1)))/2;
    0356
    0357 %figure(10); clf; plot(travdist,bin_areas,'k.-',travdist,ameth2,'r.-')
    0358
    0359 %Compute the depth for each subarea
    0360 d = bin_areas./delW;
    0361 % Fill gaps
    0362 indx1 = find(isnan(d));
    0363 indx2 = find(~isnan(d));
    0364 d(indx1) = interp1(travdist(indx2),d(indx2),travdist(indx1));  %fills the nans
    0365     % %Old depth comp
    0366     % d2=(1:sztravdist+1);
    0367     % d2(1)=ameth2(1)/midptry(1);
    0368     % d2=ameth2./(midptry-midptly);
    0369     % d2(sztravdist+1)=ameth2(sztravdist+1)./(-midptly(sztravdist+1)+travdist(sztravdist+1));
    0370     % d2 = d2';
    0371
    0372 %figure(10); clf; plot(travdist,d,'k.-')
    0373
    0374
    0375 %% Compute the Approximate Discharge
    0376 Q = nansum(lam.*bin_areas)
    0377 if Q < 0
    0378     Q = -Q;
    0379     lam = -lam;
    0380     shearvel = -shearvel;
    0381 end
    0382
    0383 %% Compute the mean cross-sectional velocity
    0384 Ua = Q./XSarea
    0385
    0386 %% Compute the LAM velocity deviation
    0387 uprime = lam - Ua;
    0388 %figure(10); clf; plot(travdist,uprime.*d,'b-')
    0389
    0390 %% Compute the transverse mixing coefficients (variable)
    0391 B = travdist(end) - travdist(1)  %Width
    0392 if 1
    0393     Ey = 0.21*shearvel.*d;  %Elder approach
    0394     % Fill gaps
    0395     indx1 = find(isnan(Ey));
    0396     indx2 = find(~isnan(Ey));
    0397     Ey(indx1) = interp1(travdist(indx2),Ey(indx2),travdist(indx1));  %fills the nans
    0398 else
    0399     % Compute the value following Shen et al. 2010
    0400     Ustar = trapz(travdist,shearvel)/B  %Cross-sectionally averaged shear velocity
    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     % Fill gaps
    0408     indx1 = find(isnan(Ey));
    0409     indx2 = find(~isnan(Ey));
    0410     Ey(indx1) = interp1(travdist(indx2),Ey(indx2),travdist(indx1));  %fills the nans
    0411 end
    0412
    0413 Eyc = trapz(travdist,Ey)/B;  %Constant value of tranverse mixing coefficient
    0414 %figure(10); clf; plot(travdist,Ey,'b-')
    0415
    0416 %% Compute the integrals (variable Ey)
    0417 % Integral 1
    0418 %indx = find(~isnan(uprime));
    0419 INT1 = cumtrapz(travdist,uprime.*d);
    0420 %INT1 = cumtrapz(travdist(:),uprime(:).*d(:));
    0421
    0422 % Integral 2
    0423 INT2 = cumtrapz(travdist,INT1./(Ey.*d));
    0424 %INT2 = cumtrapz(travdist(:),INT1./(Ey(:).*d(:)));
    0425
    0426 % Integral 3
    0427 INT3 = cumtrapz(travdist,INT2.*uprime.*d);
    0428 %INT3 = cumtrapz(travdist(indx),INT2.*uprime(indx).*d(indx));
    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 %% Compute and return the longitudinal dispersion coefficient
    0443
    0444 k = -INT3(end)/XSarea;
    0445 %percentdiff=100*(2111.7-k)/(2111.7);
    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 %disp(' ')
    0450 %disp(['Davide ended up with 2111.7 m^2/s. Which is a ' num2str(percentdiff) '% difference.  '])
    0451
    0452 %% Compute the integrals (constant Ey)
    0453
    0454 % Integral 2
    0455 INT2c = cumtrapz(travdist,INT1./(Eyc.*d));
    0456 %INT2 = cumtrapz(travdist(:),INT1./(Ey(:).*d(:)));
    0457
    0458 % Integral 3
    0459 INT3c = cumtrapz(travdist,INT2c.*uprime.*d);
    0460 %INT3 = cumtrapz(travdist(indx),INT2.*uprime(indx).*d(indx));
    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 %% Compute and return the longitudinal dispersion coefficient
    0471
    0472 kc = -INT3c(end)/XSarea;
    0473 %percentdiff=100*(2111.7-kc)/(2111.7);
    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 %% Plot a combined plot
    0479 figure(8); clf
    0480 subplot(3,1,1); plot(travdist,beddepth,'k-',travdist,d,'r-'); ylabel('Bed Depth (m)')
    0481 %subplot(4,1,2); plot(travdist_orig,zo,'k-'); ylabel('z0 (m)')
    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