
%   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
      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;
    end;


%     adjust this dependent on 60m isobath
%     find dist based 60m isobath
    if ref==5;
      load isobath60;
      disp('Hardcoded to reference 60m isobath');
    end;

%   if first lat is less than last reverse the "x" reference
    if yadcp(1,2)<yadcp(length(yadcp),2);
      cumdist=-1*(cumdist-max(cumdist)); % make distance north to south
      if ref==2;
        reflat=interp1(lon60,lat60,yadcp(length(yadcp),1));
        ddn=dist(reflat,yadcp(1,1),yadcp(length(yadcp),2),yadcp(length(yadcp),1));
        if reflat<yadcp(length(yadcp),2);ddn=-ddn;end;
      end;
    else;
      if ref==2;
        reflat=interp1(lon60,lat60,yadcp(1,1));
        ddn=dist(reflat,yadcp(1,1),yadcp(1,2),yadcp(1,1));
        if reflat<yadcp(1,2);ddn=-ddn;end;
      end;
    end;
    if ref==2;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=[.2 5]; 
    end;
    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');
    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;
       [ut,vt]=detide(yadcp(:,1),yadcp(:,2),yadcp(:,9)*24+(yr-94)*365*24+ad1);

%      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;

    zi=objmap(cumdist',ndep,yadcp(:,wv),X,Y,xysc,zto,[40 30 100 1/3 1]);

%   note: if we are detiding with quoddy we need to substract the tide
%   by first gridding the tide to the same resolution as adcp
    if dtm==2;
      zia=zi;
      disp('Detiding with QUODDY output...')

%     load modeled tides
%      load /home3/ocn/jmanning/model/hcast/bank150/ts_912may97.dat
      load /data5/jmanning/model/bank150/may97/r62/ts_912may97.dat
      ts=ts_912may97;

%     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));
      disp(['Harcoded to rotate modeled velocity ',num2str(rotate),' degrees']);
      zit=objmap(cumdistt',ndept,cr,Xt,Yt,xysc,zto,[40 30 100 1/3 1]);

%     Now we subtract the tide from the adcp 
      zi=zia-zit;
    end;      

%   convert to cm/s
    disp('Note: If you want to make west positive to show westward in blue!!!');
    zi=100*zi;
%    zi=-100*zi;
% %%%%%%%%%%%%%%%%%%  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=2;
      v=[-5:inc:20];
    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');
%    pcolor(xi,yi,zi)
%    shading interp
    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];

%   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-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;
%      kks=4;
      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 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
    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;min(wdep);min(wdep);wdep(1)]];
%      [wdep;wdep(length(wdep));wdep(1)]];
    else;
      blk=[[cumdist3';cumdist3(length(cumdist3));cumdist3(1);cumdist3(1)],...
      [wdep;min(wdep);min(wdep);wdep(1)]];
%          [wdep;dmax;dmax;wdep(1)]];
%      [wdep;wdep(length(wdep));wdep(1)]];
    end;

    pltp=get(gca,'position');
%    set(gca,'position',[pltp(1) pltp(2)+.05 pltp(3) pltp(4)]);
    if ksec==nsec;
%      caxis([floor(min(min(zi))) ceil(max(max(zi)))]);
       caxis([min(v) max(v)]);
       ci=inc;cl=inc*2;
      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;

%   add tidal phase labels
    if dtm~=4;
      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
