%  plotAvgSz.m  - plot horizontal average of salinity versus
%                 depth for a set of stations.
%                 Includes one standard deviation envelopes.
%    needs variables set:  DataDir CruiseName DepthResolution ctdcastlist 
%                          srange deprange InfoString
%    sets variables: none
%    calls: readctd

S=[10 11];
dep=[0 1];
figure;

plot(S,dep);
axis([srange deprange]);
title('Average Salinity vs depth');
xlabel('Salinity');
ylabel('Depth (m)');
text(srange(1) + .1*(srange(2)-srange(1)),.75*deprange(1),InfoString);
text(srange(1) + .1*(srange(2)-srange(1)),.81*deprange(1),CruiseName);
hold on;

%              set interpolation depths at 10 m interval
zi=deprange(1):10:deprange(2)-10;
SumS=zeros(size(zi));
SumSqS=zeros(size(zi));
NumberS=zeros(size(zi));

for cast = ctdcastlist
  file= [DataDir '/' CruiseName '/Ctd/' DepthResolution '/ctd' castnum(cast) '.dat'];
   [nlines,ndepths,desc,p,t,s,pt,r0,r1,bv,t0,t1,s0,s1,o2,par,tr,fl] = readctd(file);
  dep=-p;
  si=interp1(dep,s,zi,'linear');
  i=find(~isnan(si));
  SumS(i)=SumS(i)+si(i);
  SumSqS(i)=SumSqS(i)+si(i).^2;
  NumberS(i)=NumberS(i)+1;
end
AvgS=SumS./NumberS;
VarS=SumSqS./NumberS - (SumS./NumberS).^2;
SdS=sqrt(VarS);
lowerS=AvgS - SdS;
upperS=AvgS + SdS;

plot(AvgS,zi,'b',lowerS,zi,'r--',upperS,zi,'r--'); 

hold off;
