% 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

clf;
ca;

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

% define year
yr=97;

% define # of panels
npan=4;

subplot(npan,1,1)


% limit time window
if yr==93;
  st=141.5;et=147.5;
elseif yr==94;
  st=140.5;et=146.5;
elseif yr==95;
  sta=188;
  if sta==188;
    st=133.51;et=135;
  end;
elseif yr==97;
%  sta=82
   sta=155
  if sta==82;
    st=99.4354;
    et=105.7722;
  elseif sta==155;
    st=132.645;
    et=137.677;
  end;
end;

% buoy data
eval(['load /data5/jmanning/wind/mat/b11',num2str(yr)]);
id=find(t>=st&t<=et&airt<99);
y=y(id,:);t=t(id);
uw=y(:,1);vw=y(:,2);
% 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');
if yr==98;md=1;else;md=2;end;
if md==1;
  load /data5/jmanning/moor/moor93/m15m1hr.dat
  load /data5/jmanning/moor/moor93/m45m1hr.dat
  m1=m15m1hr;
  m2=m45m1hr;
% fit this time series to buoy
  u1=interp1(m1(:,2),m1(:,3),t);
  v1=interp1(m1(:,2),m1(:,4),t);
  mt1=interp1(m1(:,2),m1(:,7),t);
  u2=interp1(m2(:,2),m2(:,3),t);
  v2=interp1(m2(:,2),m2(:,4),t);
  mt2=interp1(m2(:,2),m2(:,7),t);
% calculate vertical average from mooring
  ua=(u1+u2)/2;
  va=(v1+v2)/2;  
  y=[airt(id)' wspd' va y];
  stkcombo(yr,t',y,1,4,2,['air temp (degC)';...
                          'wind spd (m/s) ';...
                          'N/S flow (cm/s)';...
                          'wind (m/s)     '],10)
else;
  if yr==94|yr==97;
  y=[airt(id)'; wspd; y'];
  y=y';
  elseif yr==93|yr==95;
  y=[airt(id); wspd; y'];
  y=y';
  end;
  stkcombo(yr,t',y,1,3,2,['air temp (degC) ';...
                        'wind speed (m/s)';...
                        'wind (m/s)      '],5)
end;

aaa=axis;
[axl]=get(gca,'xlim');
aap=get(gca,'position');
% make this subplot plot a little taller
aap(4)=aap(4)+.05;
set(gca,'position',aap);
aap(3)=0.64;


% 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;lat1=lat;lon1=lon;
  load /data5/jmanning/ctd/jgofout/al9306/nmfs_ctd_1993b
  lyd=[x1;lyd];press=[y1;press];temp=[z1;temp];sal=[s1;sal];cast=[c1;cast];
      lat=[lat1;lat];lon=[lon1;lon];
  yrday0_gmt=lyd+4/24-1;
elseif yr==94;
  load /data5/jmanning/ctd/jgofout/al9403/nmfs_ctd_1994all

%%%%%%%%%%%  WARNING: DOUBLE CHECK YEARDAY
  lyd=yrday1_local;

elseif yr==95;
  if sta==188;
    load /data5/jmanning/ctd/jgofout/sj9507/nmfs_ctd_1995_188
    yrday0_gmt=yrday1_gmt;
  elseif sta==197;
    load /data5/jmanning/ctd/jgofout/sj9507/nmfs_ctd_1995
  end;
elseif yr==97;
   if sta==82;
     load /data5/jmanning/ctd/jgofout/oc301/nmfs_ctd_82
   else;
      load /data5/jmanning/ctd/jgofout/oc303/nmfs_ctdsb
   end;
   temp=temp_ctd;
   sal=sal_ctd;
   clear sigmat;
end;

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

% subsample
disp([num2str(length(cast)), ' points']);
%npts=input('Number of ctd pts to skip?');
npts=2;
disp(['HARDCODED to skip every ',num2str(npts),' points']);

% correct yrday1_local to yrday0_gmt (if necessary subtract 1)
% differs each year

x=yrday0_gmt(1:npts:length(yrday0_gmt));
y=-1*press(1:npts:length(yrday0_gmt));
z=temp(1:npts:length(yrday0_gmt));
s=sal(1:npts:length(yrday0_gmt));
c=cast(1:npts:length(yrday0_gmt));
lat=lat(1:npts:length(yrday0_gmt));
lon=lon(1:npts:length(yrday0_gmt));


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

% Choose casts & build the new collection of data
%cvec=input('Enter cast number vector: ');
%if yr==93;;cvec=[6:12 40:50 53:64 159:163];
if yr==93;
  cvec=[6:12 53:64 159:163];
elseif yr==94;
  cvec=[1:11 61:73];
elseif yr==95;
  if sta==188;
    cvec=[1 3:6];
  else;
    cvec=[49:1:55 57];
  end
elseif yr==97;
  if sta==82;
    cvec=[2:14 18:19];
  else;
    cvec=[2 4:7 9:11];
  end;
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);latn=lat(id);lonn=lon(id);
  else;
    xn=[xn;x(id)];yn=[yn;y(id)];zn=[zn;z(id)];sn=[sn;s(id)];cn=[cn;c(id)];latn=[latn;lat(id)];lonn=[lonn;lon(id)];;
  end;
end;
x=xn;y=yn;z=zn;s=sn;c=cn;lat=latn;lon=lonn;

% set sc=1 to skip contours;
sc=0;

   
% 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

  if ksec==1;
    var='temperature';
    v=[5.:.2:9.];
    ci=.1;cl=.2;
  elseif ksec==2;z=sigmat(s,z)-1000;
    var='sigma-t';
    v=[25.0:.1:25.7];
    ci=.025;cl=.05;
  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,.5,[1 5]);

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

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

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


  if ptype==1;
	pcolor(yrday_gmt(1:ncasts) ,pmatrix(1:maxmax,innlonxmax),zz(1:maxmax,1:ncasts));
	hold on
	cs=contour(yrday_gmt(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;
%      wop=menu('what?','color shade','line contour');
       wop=2;
        if wop==1;
          disp('Hardcoded to color shades');
	  econtour(xi,yi,zi,v,'label',140,'fill','color','k','fontsize',10);
        else;
	  econtour(xi,yi,zi,v,'label',140,'color','g','fontsize',10);
        end;
    else;
	  econtour(yrday_gmt(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(yrday_gmt(1:ncasts),pmatrix(1:maxmax,indexmax),zz(1:maxmax,1:ncasts),v);
    end;
  end;

%%%%%%%%%%%%%%%%  Blankout bottom %%%%%%%%%%%%%%%%%%%%%

  maxy=-80;
  id=find(diff(c)~=0);
% find bottom values for this subset
  idbot=[id; length(c)];
  yblk=[y(idbot)+5; maxy; maxy; y(idbot(1))+5];
  xblk=[x(idbot); et; floor(st);floor(st)];
  hold on;
  fill(xblk,yblk,'y');


%%%%%%%%%%%%%%%%  Fix Axis & Label %%%%%%%%%%%%%%
  axis([floor(st) et maxy 0])
  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);
    id=[id; length(c)];
    skip=1;
    for i=1:skip:length(id);
      h=text(x(id(i))',3.,num2str(c(id(i))));
      set(h,'FontSize',8,'HorizontalAlignment','center')      
    end
  elseif ksec==nsec;
    ylabel('Depth(m)')
  end;
  set(gca,'xticklabels','')
  text(st+.5,-40,var);
  hgt=.17;
  if yr==93;
    set(gca,'position',[aap(1)    aap(2)-hgt*ksec+hgt    0.60    hgt])
  elseif yr==94;
    set(gca,'position',[aap(1)    aap(2)-hgt*ksec-.04    0.61    hgt])
  elseif yr==95;
    set(gca,'position',[aap(1)    aap(2)-hgt*ksec-.04    aap(3)    hgt])
  elseif yr==97;
    set(gca,'position',[aap(1)    aap(2)-hgt*ksec-.05    aap(3)    hgt])
  end;

% end of multiple contour panels
end;

% post biology on sigmat plot
  if yr<95;
    eval(['load /home3/ocn/jmanning/gbs',num2str(yr),'/gregs/strat',num2str(yr),'c.dat']);
    eval(['ss=strat',num2str(yr),'c']);
    hold on;
    coh=questdlg('What?','Cod','Haddock','both');
    dotc=questdlg('What color dot?','m','k','g');
%    wf=questdlg('Which form?','no./net','no./10m2');
    wf='no./10m2';
%   calculate gmt yrday0_gmt
    [nr,nc]=size(ss);
    sstime=julian(ss(:,3),ss(:,4),ss(:,2)-1900,0*ones(nr,1))+ss(:,5)+4/24*ones(nr,1);
    ss(:,15:17)=-1*ss(:,15:18);
    sizesym(sstime,ss(:,18),ss(:,15),2,.5,dotc);
    text(floor(st)+(et-st)/20,mean(yi),coh,'color','w','fontsize',18)
    text(floor(st)+(et-st)/30,-50,wf,'color','w','fontsize',16)

  elseif yr==95;
    eval(['load /home3/ocn/jmanning/gbs95/moc/s',num2str(sta),'.pre']);
    eval(['ss=s',num2str(sta)]);
    hold on;
    
    id=find(ss(:,4)>-100);
    coh=menu('What?','Cod','Haddock','both');
    dotc=questdlg('What color dot?','m','k','g');
    wf=questdlg('Which form?','no./net','no./100m3');
    if wf=='no./100m3';
%      standardize these
       ss(:,5)=ss(:,5)./ss(:,7)*100;
       ss(:,6)=ss(:,6)./ss(:,7)*100;
    end;
    ss(:,5:6)=-1*ss(:,5:6);
    if coh==1;
      sizesym(ss(id,1),ss(id,4),ss(id,5),2,.5,dotc);
      text(floor(st)+(et-st)/20,mean(yi),'Cod','color','w','fontsize',18)
    elseif coh==2;
      sizesym(ss(id,1),ss(id,4),ss(id,6),2,.5,dotc);
      text(floor(st)+(et-st)/20,mean(yi),'Haddock','color','w','fontsize',18)
    elseif coh==3;
      sizesym(ss(id,1),ss(id,4),ss(id,5)+ss(id,6),2,.5,dotc);
      text(floor(st)+(et-st)/20,mean(yi),'Larvae','color','w')
    end;
    text(floor(st)+(et-st)/30,-50,wf,'color','w','fontsize',16)
  end;

elseif sc==1;
% end of contouring skip
  disp('************  WARNING *************')
  disp(' HARD CODED TO NOT SKIP THE GRIDDING & CONTOURING');
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);
zz=z(idzd);
id=find(diff(c(idzd))~=0);
%   find bottom values for this subset
idbot1=[id; length(c(idzd))];
%   find surface values for this subset
idsurf1=[1;idbot1(1:length(idbot1)-1)+1];
%   subtract the bottom values from surface
for ik=1:length(idsurf1);
    diffst(ik)=zz(idbot1(ik))-zz(idsurf1(ik));
end;
%plot(x(idzd(idsurf1)),diffst,'r','linewidth',3)
%text(st+0.5,mean(diffst),'sigma-t (@20m-surface)');

% calulate max stratification parameter "sp"
id=find(diff(c)~=0);
%   find bottom values for this subset
idbot=[id; length(c)];
%   find surface values for this subset
idsurf=[1;idbot(1:length(idbot)-1)+1];
%   use surface & bottom values as reference
for ik=1:length(idsurf);
%   new modification to find only those with depths <30m
    idfc=idsurf(ik):idbot(ik);
    id30=find(y(idfc)>-30);
%   now calculate max grad in this depth range
    sp(ik)=max(diff(z(idfc(id30))));
    idzp=find(diff(z(idfc(id30)))==sp(ik));
    zp(ik)=y(idsurf(ik)+idzp);
%   now dimensionalize this correctly with the depth increment
    dinc=abs(y(idsurf(ik)+idzp)-y(idsurf(ik)+idzp+1))
    sp(ik)=sp(ik)/dinc;
end;
hold on;
%plot(x(idsurf),sp,'g--','linewidth',3)
for iz=1:length(zp);
  zpc=num2str(-1*zp(iz));
  idd=find(zpc=='.');
  text(x(idsurf(iz)),sp(iz),zpc(1:idd-1),'fontsize',6);
end;

% make a double plot with these variables
[ax,pl,ylab]=plotyy(x(idzd(idsurf1)),diffst,x(idsurf),sp,...
 '"SP" (max density grad/m <30m)',...
 [0 0 0 floor(st) et min(diffst) max(diffst) min(sp) max(sp)]);
set(pl(1),'color','r','linewidth',3);
set(pl(2),'color','g','linestyle','--','linewidth',3);
hgt=.17;
set(ax(1),'position',[aap(1)    0.55-hgt*(nsec+1)+hgt    aap(3)    hgt])
set(ax(2),'position',[aap(1)    0.55-hgt*(nsec+1)+hgt    aap(3)   hgt]);
ylabel('solid = sigma-t (@20m-surface)');

%legend('"ST" sigma-t (@20m-surface)','"SP" (max density grad/m <30m)')

%axis([floor(st) et min(sp) max(diffst)]);

%yd2calab(yr);
xlabel(['Yearday (yrday0_gmt) 19',num2str(yr)]);

orient tall

% write out parameters to a file for Greg
% (note: we want diffst, sp, zp, wind (y), wind (y(i-1:i)), taux, tauy,
%        u, v, tide at the start of each cast.)
% First, we need to get the 2 & 4 hour wind from "y" interpolated onto
%        the 2 and 4 hours before the observation times "x"
% find a distinct time, cast number, and latlon
x=unique(x);
dc=unique(c)';
idl=find(diff(c)~=0);
lat=[lat(1);lat(idl)];
lon=[lon(1);lon(idl)];
%y0_gmt=[x(1);x(idl)];
for j=1:length(x);
  for k=1:2;
    wn(j,k)=mean(interp1(t,wspd,[x(j)-(k*2)/24:1/24:x(j)]));
    un(j,k)=mean(interp1(t,uw,[x(j)-(k*2)/24:1/24:x(j)]));
    vn(j,k)=mean(interp1(t,vw,[x(j)-(k*2)/24:1/24:x(j)]));
    [taux,tauy]=wstr(un(j,k),vn(j,k),10);
    tau(j,k)=sqrt(taux^2+tauy^2);
  end;
end;

% Second, get the tide at the time of each cast or Moc tow
% com=menu('tide at what ops?','ctd','moc');
com=2;


% calculate number of years since 1994
ha94=(yr-94);
if com==1;
  [ut,vt]=detide(lon,lat,ha94*24*365+x*24);

% Now, write it all out to a file 
  yy=[dc; diffst; sp; zp; wn(:,1)'; wn(:,2)'; tau(:,1)'; tau(:,2)'; ut'; vt'; x']';
  str='       ctd     sigmat20m    maxgrad  maxgradz<30    wind2h    wind4h    wstress2h   wstress4h     u         v   yrday0_gmt' ;
  fout=['gregsn',num2str(yr),'.dat']
  fploopy(yy,fout,'a',str);
elseif com==2;
  disp('Hardcoded to output at MOCNESS times rather than CTD times');
  eval(['load /home3/ocn/jmanning/gbs',num2str(yr),'/moc/sxxx.dat']);
  mth=5;
  gmt=sxxx(:,3)+400;
  if gmt>2400;
    sxxx(:,4)=sxxx(:,4)+1;
    gmt=2400-gmt;
  end;
  ydd=julian(mth*ones(length(sxxx),1),sxxx(:,4),yr*ones(length(sxxx),1),gmt);  
% get posn from nearby ctd
  for kk=1:length(sxxx);
    idll=find(dc==sxxx(kk,2));
    mocla(kk)=lat(idll(1));
    moclo(kk)=lon(idll(1));
  end;
  [ut,vt]=detide(moclo,mocla,ha94*24*365+ydd*24);

% get wind info at moc times from buoy data
  eval(['load /data5/jmanning/wind/mat/b11',num2str(yr)]);
  id=find(t>=st&t<=et);
  y=y(id,:);t=t(id);
  uw=y(:,1);vw=y(:,2);
% generate wind speed
  clear wspd;
  for k=1:length(y);
    wspd(k)=sqrt(y(k,1)^2+y(k,2)^2);
  end;
  for j=1:length(sxxx);
%   get a two hour wind and a 4 hour wind
    for k=1:2;
      mwn(j,k)=mean(interp1(t,wspd,[ydd(j)-(k*2)/24:1/24:ydd(j)]));
      mun(j,k)=mean(interp1(t,uw,[ydd(j)-(k*2)/24:1/24:ydd(j)]));
      mvn(j,k)=mean(interp1(t,vw,[ydd(j)-(k*2)/24:1/24:ydd(j)]));
      [mtaux,mtauy]=wstr(mun(j,k),mvn(j,k),10);
      mtau(j,k)=sqrt(mtaux^2+mtauy^2);
    end;
  end;


% Now, write it all out to a file 
  yy=[sxxx mwn(:,1) mwn(:,2) mtau(:,1) mtau(:,2) ut vt ydd]
  str='       moc#     ctd#      hhmm      day   wind2h    wind4h    wstress2h   wstress4h          ut       vt     yrday0_gmt';
  fout=['/home3/ocn/jmanning/gbs',num2str(yr),'/gregsn',num2str(yr),'tide.dat']
  fploopy(yy,fout,'a',str);

end;




% attempt to calculate statsistics  
%load /home3/ocn/jmanning/gbs93/gregsn93_mat.dat
%m3=gregsn93_mat;
%mean93=mean(m3(:,1))
%load /home3/ocn/jmanning/gbs94/gregsn94_mat.dat
%m4=gregsn94_mat;
%mean94=mean(m4(:,1))
meansp=mean(sp)
meanst=mean(st)