function []=dh(cdfin)
%function cmean(cdfin)
%
% Dynamic Height overwrites "AM" variable
%    of a ecomsi model run.  cdfin is the ecomsi.cdf run, and
%    cdfout is a dummy ecomsi.cdf run that has at least three time
%    steps of varout that you don't mind writing over.  The mean value
%    is stored in the first time step, max in the 2nd time step, and
%    min in the 3rd time step 
%  
%  varin is the variable ('conc' if not supplied);
%
%   Example:  cmean('ecomsi.cdf.1','vmean.cdf.1',1:15,'salt')
%   would calculate the mean, max and min value for salt from time
%   steps 1 to 15 and store the resulting fields in vmean.cdf.1
% 
cdfid1=mexcdf('open',cdfin,'write');

[nam,nx]=mexcdf('diminq',cdfid1,'xpos');
[nam,ny]=mexcdf('diminq',cdfid1,'ypos');
[nam,nz]=mexcdf('diminq',cdfid1,'zpos');
[nam,nt]=mexcdf('diminq',cdfid1,'time');
% 
[sigma]=mexcdf('varget',cdfid1,'sigma',[0],[nz]);
dz=-diff(sigma);
% 
[depth]=mexcdf('varget',cdfid1,'depth',[0 0],[ny nx],1);
depth=replace(depth,-99999,0);

%mexcdf('REDEF', cdfid1);
%mexcdf('varrename',cdfid1,'am','dh');
%mexcdf('endef',cdfid1);

for i=1:nt,
  e=mexcdf('varget',cdfid1,'elev',[i-1  0 0],[1 ny nx],1);
% surface dh is zero
  dh=zeros(depth);
  mexcdf('varput',cdfid1,'dh',[i-1 0 0 0],[1 1 ny nx],dh,1);
  for j=[1:nz-1],
    s=mexcdf('varget',cdfid1,'salt',[i-1 j-1 0 0],[1 1 ny nx],1);
    t=mexcdf('varget',cdfid1,'temp',[i-1 j-1 0 0],[1 1 ny nx],1);
    [svan,sigt]=swstate(s,t,zeros(t));
    dh=dh+(sigt/1000+1).*depth*dz(j)/1.030;
    mexcdf('varput',cdfid1,'dh',[i-1 j 0 0],[1 1 ny nx],dh,1);
  end
end
mexcdf('close',cdfid1);
