function cthresh(cdfin,cdfout,tind,thresh)
%function cthresh(cdfin,cdfout,tind,thresh)
%
% CTHRESH computes the fraction of time during which the concentrations
%    values at each location exceed a specified threshold for a 
%    an ecomsi model run.  
%    cdfin is the ecomsi.cdf run, and
%    cdfout is a dummy ecomsi.cdf run that has at least 3
%       steps of varout that you don't mind writing over.
%  tind is the indices of time steps to average over (1 is first step)
%  thresh is a 3 element vector containing threshold values
varin='conc';
varout='conc';
cdfid1=mexcdf('open',cdfin,'nowrite');
cdfid2=mexcdf('open',cdfout,'write');

[nam,nx]=mexcdf('diminq',cdfid1,'xpos');
[nam,ny]=mexcdf('diminq',cdfid1,'ypos');
[nam,nz]=mexcdf('diminq',cdfid1,'zpos');
nn=length(thresh);
for j=[1:nz-1],
  st1=zeros(nx,ny);
  st2=zeros(nx,ny);
  st3=zeros(nx,ny);
  for i=tind,
    s=mexcdf('varget',cdfid1,varin,[i-1 j-1 0 0],[1 1 ny nx]);
    sind=s>thresh(1);
    st1=st1+sind;
    sind=s>thresh(2);
    st2=st2+sind;
    sind=s>thresh(3);
    st3=st3+sind;
  end
  st1=100*st1/length(tind);
  mexcdf('varput',cdfid2,varout,[0 j-1 0 0],[1 1 ny nx],st1);
  st2=100*st2/length(tind);
  mexcdf('varput',cdfid2,varout,[1 j-1 0 0],[1 1 ny nx],st2);
  st3=100*st3/length(tind);
  mexcdf('varput',cdfid2,varout,[2 j-1 0 0],[1 1 ny nx],st3);
  end
end
mexcdf('close',cdfid1);
mexcdf('close',cdfid2);
