function [stress,u100]=bstress(infile,cdmin,z0)
% [stress,u100]=bstress(infile,cdmin,z0)
% Calculates maximum bottom stress via the formulation in ECOM-SI using 
% velocities contained in the file infile. cdmin is the minimum drag
% coefficient.  z0 is the roughness length. This is an ECOM-SI-specific
% m-file.

if(nargin==1),
  cdmin=.0025;
  z0=.003;
end
[dummy,IM]=mcdinq(infile,'xpos');
[dummy,JM]=mcdinq(infile,'ypos');
[dummy,KB]=mcdinq(infile,'zpos');
[dummy,times]=mcdinq(infile,'time');
sigma=mcvgt(infile,'sigma',[0],[KB]);
depth=mcvgt(infile,'depth',[0 0],[JM IM]);
% NOTE: time dependence of drag coefficient is ignored in ECOM-SI
cd=max(cdmin*ones(IM,JM),.16*ones(IM,JM)./log((sigma(KB-1)-sigma(KB))/z0*depth).^2);
stress=zeros(times,IM*JM);
for i=0:times-1;
  u=mcvgt(infile,'u',[i KB-2 0 0],[1 1 JM IM]);
  u = remat(u,JM);
  ushift = zeros(u);
  ushift(1:IM-1,1:JM-1) = u(2:IM,2:JM);
  u = .5*(u+ushift);
  u = u(:)';
  v=mcvgt(infile,'v',[i KB-2 0 0],[1 1 JM IM]);
  v = remat(v,JM);
  vshift = zeros(v);
  vshift(1:IM-1,1:JM-1) = v(2:IM,2:JM);
  v = .5*(v+vshift);
  v = v(:)';
  stress(i+1,:)=cd(:)'.*(u.*u+v.*v);
end
stress=remat(max(stress),JM);
zr=1.;
cd100=z0tocd(z0,zr);   % cd at 1 m
u100=sqrt(stress/cd100);
