% routine to plot a multi panel time series of
% airtemp,sst,wind,stckwind,sigmatt,and temp 
% using the combination of shipboard, buoy,  and ctd data



path(path,'/home3/ocn/jmanning/atsee/stick')
path(path,'/home3/ocn/jmanning/atsee/saga')

% define year
yr=94;

% define # of panels
npan=4;

subplot(npan,1,1)


% limit time window
%st=141;et=147;
st=144.5;et=146.5;

% buoy data
eval(['load /data5/jmanning/wind/mat/b11',num2str(yr)]);
id=find(t>=st&t<=et);
y=y(id,:);t=t(id);

% generate wind speed
for k=1:length(y);
  wspd(k)=sqrt(y(k,1)^2+y(k,2)^2);
end;

% get mooring data?
%md=menu('Do you have mooring data?','y','n');
md=2;
if md==1;
  load /data5/jmanning/moor/moor93/m15m1hr.dat
  m=m15m1hr;
% fit this time series to buoy
  u=interp1(m(:,2),m(:,3),t);
  v=interp1(m(:,2),m(:,4),t);
  mt=interp1(m(:,2),m(:,7),t);
  stkcombo(yr,t',y,1,4,2,['air temp (degC)';...
                        'N/S flow (cm/s)';...
                        'wspd (m/s)     ';...
                        'wind (m/s)     '],10)
else;
  y=[airt(id) wspd'];
  stkcombo(yr,t',y,0,2,2,['air temp (degC)';...
                        'wspd (m/s)     ']);
end;

aaa=axis;
aap=get(gca,'position');



% sometimes need two files
if yr==93;
  load /data5/jmanning/ctd/jgofout/al9306/nmfs_ctd_1993a
  x1=lyd;y1=press;z1=temp;s1=sal;c1=cast;
  load /data5/jmanning/ctd/jgofout/al9306/nmfs_ctd_1993b
  lyd=[x1;lyd];press=[y1;press];temp=[z1;temp];sal=[s1;sal];cast=[c1;cast];
else;
  load /data5/jmanning/ctd/jgofout/al9403/nmfs_ctd_1994
end;

%%%%%%%%%%%%%%%%%%   Delimit CTD Data  %%%%%%%%%%%%%%%%%%

% subsample
disp([num2str(length(lyd)), ' points']);
npts=input('Number of ctd pts to skip?');
% correct yrday1_local to yrday0_local
x=lyd(1:npts:length(lyd))-1;
y=-1*press(1:npts:length(lyd));
z=temp(1:npts:length(lyd));
s=sal(1:npts:length(lyd));
c=cast(1:npts:length(lyd));

% rid of NaN's
id=find(isnan(x)~=1);
x=x(id);y=y(id);z=z(id);c=c(id);s=s(id);

% Choose casts & build the new collection of data
%cvec=input('Enter cast number vector: ');
%if yr==93;;cvec=[6:12 53:64 159:163];elseif yr==94;cvec=[1:7 61:73 79:82];end
%if yr==93;cvec=[ 53:64 ];elseif yr==94;cvec=[1:7 61:73];end
if yr==93;cvec=[ 53:64 ];elseif yr==94;cvec=[61:73];end
for k=1:length(cvec);
  id=find(c==cvec(k));
  if k==1;
    xn=x(id);yn=y(id);zn=z(id);sn=s(id);cn=c(id);
  else;
    xn=[xn;x(id)];yn=[yn;y(id)];zn=[zn;z(id)];sn=[sn;s(id)];cn=[cn;c(id)];
  end;
end;
x=xn;y=yn;z=zn;s=sn;c=cn;
   
% Now the nsec contours
nsec=2;  
for ksec=1:nsec;
  subplot(npan,1,ksec+1);

% Change parameters relative to section #
% make it density in second section
% choose contour levels and cbar levels
% ci=(max(z)-min(z))/5;
% cl=ci;
% v=[min(z):(max(z)-min(z))/10:max(z)];

  if ksec==1;
    var='temperature';
    v=[6:.25:8.5];
    ci=.25;cl=.5;
  elseif ksec==2;z=sigmat(s,z)-1000;
    var='sigma-t';
    v=[25:.1:26];
    ci=.1;cl=.2;
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   GRID DATA  &&&&&&&&&&&&&&&&&&&&&&&

% gd=menu('Grid data?','use griddata',...
%                         'use sagas mincurvi',...
%                         'use sagas kriging',...
%                         'use sagas objmap',...
%                          'none');       
  gd=4;
 
  
  [zi,xi,yi]=gridopts(x,y,z,gd);

%%%%%%%%%%%%%%%%%%%%% DEFINE METHOD OF CONTOURING %%%%%%%%%%%%%%%%%%

%    ptype=menu('Plot type?','pcolor & contour','econtour','contourf');
  ptype=2;

%%%%%%%%%%%%%%%%%%%%%%%%%%  Do the Contouring Now %%%%%%%%%%%%%%%%


  if ptype==1;
	pcolor(lyd(1:ncasts) ,pmatrix(1:maxmax,innlonxmax),zz(1:maxmax,1:ncasts));
	hold on
	cs=contour(lyd(1:ncasts) ,pmatrix(1:maxmax,innlonxmax),zz(1:maxmax,1:ncasts),v);
	clabel(cs,'manual')
    shading interp
  elseif ptype==2;
    disp('Hardconlond to use EXTCONTOUR.');
    if gd<=4;
	  econtour(xi,yi,zi,v,'fill');
%	  econtour(xi,yi,zi,v,'label');
    else;
	  econtour(lyd(1:ncasts),pmatrix(1:maxmax,innlonxmax),zz(1:maxmax,1:ncasts),v,'fill');
    end;
%    caxis([min(v) max(v)]);
%    cbar(gca,'h',ci,cl);
  elseif ptype==3;
    disp('Hardconlond to use CONTOURF.');
    cbarjm(gca,'v',ci,cl,v);
    if gd<=4;
	  contourf(xi,yi,zi,v);
    else;
%         note that contourf does not accept NaNs so we set bottm temp to
      for kl=1:ncasts;
        idl=find(isnan(zz(:,kl))==1);
        zz(idl,kl)=zz(maxpindex(kl),kl)*ones(length(idl),1);
      end;
      contourf(lyd(1:ncasts),pmatrix(1:maxmax,indexmax),zz(1:maxmax,1:ncasts),v);
    end;
  end;


%%%%%%%%%%%%%%%%  Fix Axis & Label %%%%%%%%%%%%%%
  axis([st et -80 0])
  hgt=.17;
  set(gca,'position',[0.24    0.55-hgt*ksec+hgt    0.55    hgt])
  if ksec==1;
    h=text(mean(x),10,'CTD Stations','horizontalalignment','center');
    disp('Hardcoded to title only top most section')
    id=find(diff(c)~=0);
    skip=2;
    for i=1:skip:length(id);
      h=text(x(id(i))',2.,num2str(c(id(i))));
      set(h,'FontSize',6,'HorizontalAlignment','center')      
    end
  elseif ksec==nsec;
    ylabel('Depth(m)')
  end;
  set(gca,'xticklabels','')
  text(st+.5,-40,var);
end;

% finally sigma-t diff plot
subplot(npan,1,1+nsec+1);
%   calculate sigmat diff from surface to zd
zd=-20;
idzd=find(y>=zd);
%   find bottom values for this subset
idbot=[id; length(c)];
%   find surface values for this subset
idsurf=[1;idbot(1:length(idbot)-1)+1];
%   subtract the bottom values from surface
for ik=1:length(idsurf);
    diffst(ik)=z(idbot(ik))-z(idsurf(ik));
end;
plot(x(idsurf),diffst,'r','linewidth',3)
text(st+0.5,mean(diffst),'sigma-t (@20m-surface)');
axis([st et min(diffst) max(diffst)]);
yd2calab(yr);
set(gca,'position',[0.24    0.55-hgt*(nsec+1)+hgt    0.55    hgt])
orient tall


  
