
%   here's where we grid & contour adcp sections
%   after having run "versec.m" to setup page etc.
    if field==7|field=='East';wv=1;cfield='East';elseif field==8; wv=2;cfield='North';end;

    ref=1;
    if ref==1; % relative to first station (north to south)
      cumdist=distcalc(yadcp(:,2),yadcp(:,1));
    elseif ref==3;
%      rla=input('Enter Reference Lat (DDMM.MM):');
%      rlo=input('Enter Reference Lon (positive DDMM.MM):')
      rla=4109.00;
      rlo=6746.00;
      disp('Hardcoded to reference particular');
      [rla,rlo]=ddmm2dec(rla,rlo);
      for i=1:length(yadcp);
        cumdist(i)=dist(rla,rlo,yadcp(i,2),yadcp(i,1));
      end;
    elseif ref==5;
      %     find dist of northern most station from 60m isobath
      disp('Hardcoded to reference 60m isobath');
      load isobath60;
      disp('Hardcoded to reference 60m isobath');
      if yadcp(1,2)<yadcp(length(yadcp),2);
        reflat=interp1(lon60,lat60,yadcp(length(yadcp),1));
        ddn=dist(reflat,yadcp(length(yadcp),1),yadcp(length(yadcp),2),yadcp(length(yadcp),1));
      else;
        reflat=interp1(lon60,lat60,yadcp(1,1));
        ddn=dist(reflat,yadcp(1,1),yadcp(1,2),yadcp(1,1));
      end;
      cumdist=distcalc(yadcp(:,2),yadcp(:,1));
    end;

%   if first lat is less than last reverse the "x" reference 
%   and if reflat< northmost lat change the sign of ddn, the dist to 60m
    if yadcp(1,2)<yadcp(length(yadcp),2);
      cumdist=-1*(cumdist-max(cumdist)); % make distance north to south
      if ref==5&reflat<yadcp(length(yadcp),2);ddn=-ddn;end;
    else;
      if ref==5&reflat<yadcp(1,2);ddn=-ddn;end;
    end;

    if ref==5;cumdist=cumdist+ddn;end;

%   define depth
    ndep=yadcp(:,3);
    if wctd<10;
%      ndep=-1*ndep; % make depth negative
    end

%   set the following to the same as "nsec" in quoddy method
    ngp=25;
    disp(['Hardcoded to ',num2str(ngp),' by ',num2str(ngp),' grid.']); 
    if vxmax==1;
      xmax=input('Maximum distance (km) alongtack? ');
      xmin=input('Minimum distance (km) alongtack? ');
    end;
    if vdmax==1;
%     define maximum depth and xdist desired
      dmax=input('Enter Maximum Depth: ');
      if dmax>0;dmax=-dmax;end;
    end;
    x=1;
    if x==0;
      cumdist=yadcp(:,9);
      disp('Hardcoded to use time as x-axis');
%     Here's where we set the gridding parameters
%     zto=input('Enter some number betw 0 and 1: ');
%     xysc=input('Enter vector w/x&y scales: ');
      zto=.5;xysc=[.005 10];
    else;
      disp(['Hardcode to xmax=',num2str(xmax),' dmax=',num2str(dmax)]);
      zto=.5;xysc=[.1 5]; 
%      zto=.5;xysc=[5 5]; 
    end;
    disp(['x&y scales: ',num2str(xysc(1)),' and ',num2str(xysc(1))]);
    xi=[min(cumdist):(max(cumdist)-min(cumdist))/ngp:max(cumdist)];
    yi=[minnnan(ndep):(maxnnan(ndep)-minnnan(ndep))/ngp:maxnnan(ndep)]';
    [X,Y]=meshgrid(xi,yi);


%   note: we want to grid a certain column variable depending on the    detiding
%   note: wv+6 is the filtered velocity in predone julios case

%    dtm=menu('Detiding method?','already done',
%     'quoddy','julios','none','Moody');
    dtm=4;
    if dtm==1;wv=wv+6;elseif dtm>=2;wv=wv+4;end;

    if dtm==3;
       disp('Hardcoded to detide with Julios routine')
%      Julio's detiding done here
%      account for leapyear
       if yr>96;ad1=24;else;ad1=0;end;    % corrects for yearday
       if yr==98;sb1=1/24;else;sb1=0;end; % corrects for apparent phase problem
       [ut,vt]=detide(yadcp(:,1),yadcp(:,2),yadcp(:,9)*24+(yr-94)*365*24+ad1+sb1);

%      determine the phase of the tide for labeling
       [tideph1]=gettidph(ut(1),vt(1));
       [tideph2]=gettidph(ut(length(ut)),vt(length(vt)));
       yadcp(:,5)=yadcp(:,5)-ut;
       yadcp(:,6)=yadcp(:,6)-vt;  
    end;

    if dtm==5;
       disp('Hardcoded to detide with Moody')
%      Moody's detiding done here
       [ut,vt,ydm]=moody(3,[1 2 3 4 5],[98 5 5 0200],72);

%      determine the phase of the tide for labeling
       [tideph1]=gettidph(ut(1),vt(1));
       [tideph2]=gettidph(ut(length(ut)),vt(length(vt)));

%      interpoloate moody tide to yadcp time
       disp('Interpolating moody tide to yadcp time');
       uty=interp1(ydm,ut,yadcp(:,9));
       vty=interp1(ydm,vt,yadcp(:,9));
       yadcp(:,5)=yadcp(:,5)-uty;
       yadcp(:,6)=yadcp(:,6)-vty;  
    end;
    
    disp('gridding observed alongbank adcp');
    zi=objmap(cumdist',ndep,yadcp(:,wv),X,Y,xysc,zto,[40 30 100 1/3 1]);
    zia=zi;
     if wctd~=12;
      disp('gridding observed crossbank adcp after negative factor applied');
      ziv=objmap(cumdist',ndep,-1*yadcp(:,wv+1),X,Y,xysc,zto,[40 30 100 1/3 1]);
      disp('gridding observed vertical adcp (assumes col 10 for w)');
      ziw=objmap(cumdist',ndep,yadcp(:,wv+3),X,Y,xysc,zto,[40 30 100 1/3 1]);
      ziav=ziv;
      ziaw=ziw;
    end;

%   note: if we are detiding with quoddy we need to subtract the tide
%   by first gridding the tide to the same resolution as adcp


    if dtm==2;

      disp('Detiding with QUODDY output...')

%     load modeled tides
      load /data5/jmanning/model/bank150/may97/r62/ts_912may97.dat
      ts=ts_912may97;
%      load /data5/jmanning/model/bank150/may98/r71/ts.dat

%     extract lines for this transects time range
      idt=find(ts(:,6)>=miny&ts(:,6)<=maxy);
      ts=ts(idt,:);
      ts=ts(1:2:length(ts),:);

%     we want to process this exactly as we did the adcp data
      if ref==1|ref==5; % relative to first station
        cumdistt=distcalc(ts(:,8),ts(:,7));
      elseif ref==3;
        disp('Hardcoded to reference particular');
        [rla,rlo]=ddmm2dec(rla,rlo);
        cumdistt(i)=dist(rla*ones(length(ts),1),rlo*ones(length(ts),1),...
         ts(:,8),ts(:,7));
      end;

%     if first lat is less than last reverse the "x" reference
      if ts(1,8)<ts(length(ts),8);
        cumdistt=-1*(cumdistt-max(cumdistt)); % make distance north to south
        [tideph2]=gettidph(ts(1,13),ts(1,14));
        [tideph1]=gettidph(ts(length(ts),13),ts(length(ts),14));
      else;
        [tideph1]=gettidph(ts(1,13),ts(1,14));
        [tideph2]=gettidph(ts(length(ts),13),ts(length(ts),14));
      end;
      [tidephm]=gettidph(ts(floor(length(ts)/2),13),ts(floor(length(ts)/2),14));
      
      if ref==5;cumdistt=cumdistt+ddn;end;
     
      ndept=ts(:,12); % make depth negative
      xit=[min(cumdistt):(max(cumdistt)-min(cumdistt))/ngp:max(cumdistt)];
      yit=[minnnan(ndept):(maxnnan(ndept)-minnnan(ndept))/ngp:maxnnan(ndept)]';
      [Xt,Yt]=meshgrid(xit,yit);
%     here's where we rotate the modeled u&v to cross&along?
      [cr,al]=rot(rotate,ts(:,13),ts(:,14));
%      cr=ts(:,13);al=ts(:,14));
      disp(['Harcoded to rotate modeled velocity ',num2str(rotate),' degrees']);
      disp('gridding alongbank modeled tide');
      zit=objmap(cumdistt',ndept,cr,Xt,Yt,xysc,zto,[40 30 100 1/3 1]);
      disp('gridding crossbank modeled tide after negative factor applied');
      zitv=objmap(cumdistt',ndept,-1*al,Xt,Yt,xysc,zto,[40 30 100 1/3 1]);

%     Now we subtract the tide from the adcp
      zi=zia-zit;
      ziv=ziav-zitv;

    end;      

    % now compare before & after
    chck=0;
    if chck==1;
      figure(2);
      if dtm~=4; 
        subplot(3,1,1)
      end;
      cc=contour(xi,yi,zia);
      clabel(cc);
      title('raw');
      for i=1:length(cumdist);
        text(cumdist(i),ndep(i),num2str(yadcp(i,wv),1),'fontsize',6);
      end;
      axis([min(cumdist) max(cumdist) -80 0])
       title(['Section # ',num2str(ksec)]);
      if dtm<4;
        subplot(3,1,2)
        cc=contour(xi,yi,zit);
        clabel(cc);
        title('tide');
        for i=1:length(cumdistt);
          text(cumdistt(i),ndept(i),num2str(cr(i),1),'fontsize',6);
        end;
        axis([min(cumdist) max(cumdist) -80 0])
      end;
      if dtm~=4;
        subplot(3,1,3)
        cc=contour(xi,yi,zi);
        clabel(cc);
        title('detided');
        for i=1:2:length(xi);
          for j=1:2:length(yi)
            text(xi(i),yi(j),num2str(zi(i,j),1),'fontsize',4);
          end;
        end;
        axis([min(cumdist) max(cumdist) -80 0])
      end;
      pt=menu('print this?','y','n');
      if pt==1;
        keyboard;
        pthis;
      end;
      figure(1);
    end;

    if wctd~=12;
%   convert to cm/s & exaggerate the vertical by 1000
      zi=100*zi;
      ziv=100*ziv;
      disp('Hardcoded to NOT exaggerate the vertical velocity by 1000X');
      ziw=100*ziw;
    end;

% %%%%%%%%%%%%%%%%%%  DEFINE   CONTOUR INTERVAL %%%%%%%%%%%%%%%%%%%

    if ksec==1|vincr==1;
      incr=menu('What contour interval?',...
               'min to max by range/5',...
               'prompt user',...
               'default for particular variable',...
               'same as one but integerized',...
               'from "stan.dat"');
    end;
%    [rg,mi,ma]=rangejm(min(zi));
    disp(['range: ',num2str(min(min(zi))),' to ',num2str(max(max(zi)))]);
    if incr==1;
      v=[mi:rg/5:ma];
    elseif incr==2;
      v=input('Enter min:inc:max contour int: ');
      inc=mean(diff(v));
    elseif incr==3;
      disp('Hardcoded contour intervals to user defined to:');
      inc=4;
      v=[-30:inc:10];
    elseif incr==4;
%     attempt to intgerize the intervals
      if rg>2;
        v=floor(minnnan(zi)):floor(rg/5):ceil(maxnnan(zi));
%        eval(['v=floor(min',cfield,'):floor(mean(diff(:ceil(max',cfield,');']);
      else;
        v=[mi:rg/5:ma];
      end;
    elseif incr==5;
      load stan.dat
      v=[stan(ksec,3):stan(ksec,4):stan(ksec,5)];
    end;
    ci=mean(diff(v));cl=ci*2;
    disp(['hardcoded contour intervals at ',num2str(v)]);


    econtour(xi,yi,zi,v,'fill');
    colormap(1-jet)
    hold on;


%   restrict axis to water column depth
%   (note: the "water_depth" lines were commented out for bb2mat option
%    and substituted with "bd" 
    bd=-1*yadcp(:,4);
%    bd=-1*water_depth;     

%   restrict axis to data
    a=[min(cumdist) max(cumdist) max(bd) dmax];

% vecplot cross-bank
    xx=reshape((xi'*ones(1,26))',26*26,1);
    yy=reshape((yi*ones(1,26)),26*26,1);
    zzx=reshape(ziv,26*26,1);
    zzy=reshape(ziw,26*26,1);
    pause(1);
    disp('WARNING: cross-bank vector plot  commented out in contadcp.m')
%    skp=10;
%    vecplot2(xx(1:skp:length(xx)),yy(1:skp:length(xx)),...
%      zzx(1:skp:length(xx)),zzy(1:skp:length(xx)),100,'cm/s',2,-2,'k');
    
%    keyboard;

%   blank out data below max depth threshold
    idmd=find(diff(yadcp(:,3))>0);
%   make maximum depth penetration of adcp data
%   smooth this max depth vector some?
    disp('Smoothing Max Depth of Penetration...')
    [smd,ff]=runave(yadcp(idmd,9),ndep(idmd),15,0,0,0);

%   add a search radius worth of extrapolated data in the vertical
    smd=smd-0*xysc(2); 

    blk2=[[cumdist(idmd)';cumdist(idmd(length(idmd)));cumdist(idmd(1));cumdist(idmd(1))],...
    [smd;a(4);a(4);smd(1)]];

%     blank max depth penetration
    hold on;
%    ab=fill(blk2(:,1),blk2(:,2),'k');
%    set(ab,'edgecolor','none');
    
    if ksec==1;
      over=menu('load another data and overlay?','y','n');
    end;
    if over==1;
      kks=ksec;
      wsec=input('What section is this?');
      kks=wsec;
      eval(['load sigmam',num2str(kks)]);
%     contourlabel fonsize dependent on number of sections
      clfs=15-nsec;
      econtour(xxi,yyi,zzi,[25.2:.1:25.8],'label','linewidth',2,'color','g','fontsize',clfs);

%     now add the CTD locations & bongo haul catch as saved in first pass
%      eval(['load sigmas',num2str(kks)]);
%          for i=1:ncasts
%            h=text(cumdist2(i)',0.01*maxp(indexmax),num2str(snum(i)));
%            h=text(cumdist2(i)',0,num2str(snum(i)));
%            set(h,'FontSize',8,'HorizontalAlignment','center',...
%                 'verticalalignment','bottom')     
%          end
%          posthyd(snum,cumdist2); 
    end;

%     get bottom depth from Bob' routine OR data
%     disp('...getting bottom depth from Bobs routine')
%   wdep=-1*w_depth(yadcp(:,1),yadcp(:,2));

%    wdep=-1*water_depth;
    wdep=bd;

%    cumdist3=distcalc(lat,lon);      
    cumdist3=cumdist;

%   make bottom topo blanking
%    blk=[[cumdist3';cumdist3(length(cumdist3));cumdist3(1);cumdist3(1)],...
%    [wdep;a(4);a(4);wdep(1)]];
    idb=find(isnan(wdep)~=1);
    cumdist3=cumdist3(idb);
    wdep=wdep(idb);
    if yadcp(1,2)<yadcp(length(yadcp),2);
%      blk=[[cumdist3';cumdist3(length(wdep));cumdist3(length(wdep))],...
      blk=[[cumdist3';cumdist3(length(wdep));cumdist3(1);cumdist3(1)],...
      [wdep;dmax;dmax;wdep(1)]];
%      [wdep;wdep(length(wdep));wdep(1)]];
    else;
      blk=[[cumdist3';cumdist3(length(cumdist3));cumdist3(1);cumdist3(1)],...
      [wdep;dmax;dmax;wdep(1)]];
    end;

    pltp=get(gca,'position');
%    set(gca,'position',[pltp(1) pltp(2)+.05 pltp(3) pltp(4)]);
    if ksec==nsec|vfield==1|vincr==1;
%      caxis([floor(min(min(zi))) ceil(max(max(zi)))]);
       caxis([min(v) max(v)]);
       ci=inc;cl=inc*2;
      if vfield==1|vincr==1;
        hbar=cbar(gca,'v',ci,cl);
      else;
        hbar=cbar(gca,'h',ci,cl);
        ha=get(hbar,'position');
        set(hbar,'position',[ha(1) ha(2)-ha(2)/2 ha(3) ha(4)]);
      end;
    end;

%   add tidal phase labels
    alab=0;
    if alab==0;disp('Note: tide phase labels are commented out.');end;
    if dtm~=4&alab~=0;
      text(a(1),0,tideph1,'horizontalalignment','left',...
       'verticalalignment','top','fontsize',12);
      text(a(2)/2,0,tidephm,'horizontalalignment','center',...
       'verticalalignment','top','fontsize',12);
      text(a(2),0,tideph2,'horizontalalignment','right',...
       'verticalalignment','top','fontsize',12);
    end;

%   call subroutine to post drifters on vertical sections
%    if ksec>2;
%      postdrft;
%    end;
    disp('Hardcoded to NOT run postdrft or postmoc');
     postmoc

