%
% VERSEC.M is an attempt to contour oceanographic cross-sections
% with several options on:
%
%            input data format
%            gridding method
%            contouring method

%***********     ESTABLISH SOME DEFAULT PARAMETERS BEFORE LOOP *********

disp('Warning: This version of "versec.m" has several hardcodes to minimize prompting.')
% set data path
pth='';
% mm=input('Which option?','animation','subplots for hardcopy','fullplot');
mm=3;
%  open list.dat file in the case of movie making
if mm==1; disp('Hardcoded to make movie.');fid2=fopen('list.dat','w');
  elseif mm==2;disp('Hardcoded to make subplots.')
end;
%yr=input('2-dgit year?');
yr=99;
disp(['hardcoded to Year=',num2str(yr)]);
%nsec=input('How Many Sections? ');
nsec=1;
disp(['Hardcoded to plot ',num2str(nsec),' sections']);
% set bottom depth default values
bdop=2;%bottom depth source 1=in file, 2= add 1 meter, 6=beardsley routine
ref=1; %reference for section 1= 1st cast,4=in file; 5=60m,
wctd=7;%data type 4=data, 6=nrv, 2=ascw/outhead, 7=911w/head
field=1;%variable
cfield='Temp';%variable string
incr=2;%2=prompt; 3=hardcode; contour interval specification option
ptype=2;%1=pcolor, 2=fill, 3=contour lines
sgf=0;%save gridded field? set to 1 if yes
xmax=12;
xmin=0;
dmax=-70;
sw1=45.;  % starting shelfwidth
shlfwdth=sw1;
%
% Recent Modifications:
%
% 22 Nov 1996  - add SAGA option to gridding
% 30 Jan 1997  - added SEABIRD 911 cnv file option for input
%  4 Feb 1997  - added option to contour adcp data
%  6 Mar 1997  - added option to annimate "quoddy" (modeled) sections
% 12 Mar 1997  - added Charlie Flagg CODAS option for input
%  8 Sep 1997  - added detiding of adcp w/ QUODDY2VRVT output
% 30 Oct 1997  - added option to contour bb2mat adcp output
% 28 Jan 1998  - added "postdrft" call to post drifters & scale x position
% 
%************************************************
%
% Note #1: An alternative to this routine is to use Dick Limeburners's
% "VSEC.M" which has the advantage of reading directly the .prs header for lat/lon
% info and the disadvantage of assuming that info is correct.)
%
% Note #2: You will notice that many of the user prompts have been
% commented out and hardcoded so I do not have to answer a lot of
% questions. You may want to either uncomment them or hardcode
% the option to your liking.
%
% Note #3: There are several phases to this routine and some of them
% are nearly 100 lines each. Each phase is denoted by a ****%%%%%
% comment line. There is one large loop to do multiple sections and
% for each section there is the following phases:
%
%    -load station numbers and geographic positions
%    -calculate alongtrack distances
%    -load ctd data from individual files
%    -fill the variable matrix for pressure, temp, salt, etc.
%    -create the blanking polygon (to blank out the bottom)
%    -grid
%    -contour
%    -fix up axis and text
%
%


%%%%%%%%%%%%%%%%%%%%% LOOP FOR MULITPLE SECTIONS %%%%%%%%%%%%%%%%%%%%%
%

for ksec=1:nsec;

  if mm==1;
    clf;
    orient portrait
    vbdop=0;vref=0;vfield=0;vwctd=0;vincr=0;vsgf=0;vxmax=0;vdmax=0;vptype=0;
  elseif mm==2;
    subplot(nsec,1,ksec);
  end;

%%%%%%%%%%%%%%%%%%%%%  USER SPECIFYS OPTIONS BY CLICKING BUTTONS %%%%%

% skip by user prompts in case of movie making
  if ksec==1|mm>=2;

%   initialize the codes for changes or no changes
  
%   Add some buttons to figure one to allow user options
    bbdop=uicontrol('style','checkbox','string',['bot depth=',num2str(bdop)],'Position',   [ 10 20 150 20]);
    bref=uicontrol('style','checkbox','string',['x-reference=',num2str(ref)],'Position',   [160 20 150 20]);
    bfield=uicontrol('style','checkbox','string',['variable=',cfield],'Position',    [300 20 150 20]);
    bwctd=uicontrol('style','checkbox','string',['input data=',num2str(wctd)],'Position',[ 10 50 150 20]);
    bincr=uicontrol('style','checkbox','string',['cont int option=',num2str(incr)],'Position',  [160 50 150 20]);
    bsgf=uicontrol('style','checkbox','string',['save grid=',num2str(sgf)],'Position',[300 50 150 20]);
    bxmax=uicontrol('style','checkbox','string',['xmax=',num2str(xmax)],'Position',[10 80 150 20]);
    bdmax=uicontrol('style','checkbox','string',['dmax=',num2str(dmax)],'Position',[160 80 150 20]);
    bptype=uicontrol('style','checkbox','string',['ptype=',num2str(ptype)],'Position',[300 80 150 20]);

%
    disp(' ')
    disp('')
    disp(' ***************************************************')
    disp('     Click  the variables you want to work with. ')
    disp(' ')
    disp('     Then click and press return in this window  when done. ')
    disp(' ***************************************************')
    disp(' ')
    disp('')
    pause;
    disp('loading data ...');

    set(bbdop,'Visible','off');
    set(bref,'Visible','off');
    set(bfield,'Visible','off');
    set(bwctd,'Visible','off');
    set(bincr,'Visible','off');
    set(bsgf,'Visible','off');
    set(bxmax,'Visible','off');
    set(bdmax,'Visible','off');
    set(bptype,'Visible','off');

    vbdop=get(bbdop,'Value');
    vref=get(bref,'value');
    vfield=get(bfield,'Value');
    vwctd=get(bwctd,'value');
    vincr=get(bincr,'value');
    vsgf=get(bsgf,'Value');
    vxmax=get(bxmax,'Value');
    vdmax=get(bdmax,'Value');
    vptype=get(bptype,'Value');
  end;

%%%%%%%%%%%%%%%%%%%%%  LOAD STATION DATA %%%%%%%%%%%%%%%%%%%%%%%%%


% Here we have a different ways of loading data depending
% on the format the input data is stored. There are multiple options.

  if vwctd==1;
  clear posn;
  wctd=menu('Which data format?',...
  '1) General Oceanics',...
  '2) model 19 Seabird',...
  '3) Old NODC "report" output',...
  '4) OAX "data" output',...
  '5) JGOFS .jgf file',...
  '6) Dartmouth Model "nrv or rvt"',...
  '7) model 911 Seabird',...
  '8) Candela adcp',...
  '9) Flagg adcp',...
  '10) bb2mat adcp',...
  '11) JGOFS .mat file',...
  '12) Wiebe "curtain" file');
  end;

% setup filename reference for wctd=6 option
  if wctd==6&ksec==1;fd=7;fh=0;end;

%%%%%%%%%%%%%%%%%%%%%%%%%  GET CTD DATA  %%%%%%%%%%%%%%%%%%%%%

% the objective of the following routine 
% is to get position data and station numbers (snum,q,ncast,etc.)
  getctd;

% if this is adcp note that we skip many phases starting now
  if wctd~=9&wctd~=10&wctd~=12;

  %%%%%%%%%%%%%%%%%%%%%%%%%  ALONGTRACK DISTANCES %%%%%%%%%%%%%%%%%%%%%%%

  % determine along-track distances in kms
  % by calling a mfile called "dist"
    if vref==1;
      ref=menu('Reference Point?',...
                    'first ctd',...
                     'last ctd',...
           'particular lat/lon',...
          'distance in file itself',...
          'relative to 60m isobath');
      if ref==3;
        rlo=input('Reference longitude?');
        rla=input('Reference latitude?');
      end;
    end;     
%
%
  if ref==1;
    for i=1:ncasts;
      cumdist(i)=dist(q(1,1),q(1,2),q(i,1),q(i,2));
    end
  elseif ref==2;
    for i=1:ncasts;
      cumdist(i)=dist(q(ncasts,1),q(ncasts,2),q(i,1),q(i,2));
    end
  elseif ref==3;
%   rla=input('Enter Reference Lat (DDMM.MM):');
   rlo=input('Enter Reference Lon (positive DDMM.MM):')
    disp('Hardcoded to reference particular');
    rla=4106.00;
%    rlo=6800.00;
%    rla=4100.6;
%    rlo=6734.8;
    disp(['Hardcoded ref position: ',num2str(rla),' ',num2str(rlo)]);
    [rla,rlo]=ddmm2dec(rla,rlo);
    for i=1:ncasts;
      cumdist(i)=dist(rla,-1*rlo,q(i,1),-1*rlo);
      if q(i,1)>rla;
        cumdist(i)=cumdist(i)*(-1);
      end;
    end

  elseif ref==4;
    disp('Hardcoded to get alongtrack reference from data');
    cumdist=posn(id,2)/1000;
%    cumdist=-1*(cumdist-max(cumdist));
  elseif ref==5;
    
    for i=1:ncasts;
      cumdist(i)=dist(q(1,1),q(1,2),q(i,1),q(i,2));
    end
    disp('Hardcoded to compute x-distance from the 60m isobath');
    load isobath60;
    reflat=interp1(lon60,lat60,q(1,2));
    dd=dist(reflat,q(1,2),q(1,1),q(1,2));
    if reflat<q(1,1);dd=-dd;end;
    cumdist=cumdist+dd;
  end;   


%%%%%%%%%%%%%%%%%%%%%%%%%% LOAD DATA %%%%%%%%%%%%%%%%%%%%%%%%%%


% note we skip this section for ctd option 4
  if wctd<4;

%   load CTD files
    %ty='.pre';
     ty='.edt';
    if wctd==1;
      id='oc22d';
    elseif wctd==2;
      id='sbe';
    elseif wctd==3;
      id='al8209s';
    end;

  % Next we generate cast number & open files
  % so that we may have decreasing cast numbers if needed

    disp('Hardcoded without leading zeros in cast#')
    disp('Hardcoded for the dos forward slash')
    disp('Hardcoded to expect p,t,s,Ox,fl,tr in that order')
    for i=1:ncasts
      ic=begin+i*inc+(-1)*inc;
      if ic<=9
        eval(['load ',pth,'\',id,'00',num2str(ic),ty])
      elseif begin+i-1>9&ic<100;

        eval(['load ',pth,'\',id,'0',num2str(ic),ty])
  %     remove the zeros in the case of wctd==3;
  %     eval(['load ',pth,'\',id,num2str(ic),ty])
      else
        eval(['load ',pth,'\',id,num2str(ic),ty])
      end
    end
  

  % determine number of rows per CTD
    len=[];
    for i=1:ncasts
      ic=begin+i*inc+(-1)*inc;
      if ic<=9
  	eval(['s=size(',id,'00',num2str(ic),');']);
  	len=[len;s(1,1)];
      elseif ic>9&ic<100;
  	eval(['s=size(',id,'0',num2str(ic),');']);
  	len=[len;s(1,1)];
      else
	eval(['s=size(',id,num2str(ic),');']);
	len=[len;s(1,1)];
      end
    end

  elseif wctd==4;
  % determine number of rows per CTD
    len=[];
    for i=1:ncasts;
      ic=snum(i);
      idd=find(posn(:,11)==ic&posn(:,3)<=dmax*(-1));
      len=[len;length(idd)];
    end;

  elseif wctd==5;
    id='sj07d';
    ty='.dat';

  % Next we generate cast number & open files
  % so that we may have decreasing cast numbers if needed

    disp('Hardcoded without leading zeros in cast#')
    disp('Hardcoded for the UNIX forward slash')
    disp('Hardcoded to expect p,t,s,Ox,fl,tr in that order')
    len=zeros(ncasts,1);
    for i=1:ncasts
      ic=begin+i*inc+(-1)*inc;
      if ic<=9,zer='00';elseif ic>9&ic<100;zer='0';else;zer='';end;
      eval(['idf=fopen(''',pth,'/',id,zer,num2str(ic),ty,''');'])

  % determine number of rows per CTD & get the data
      
      for k=1:9;
        lll=fgetl(idf);
      end;
      res=0;
      dd=[];
      while res~=-1&res~=[];
        res=fscanf(idf,'%f %s\n',13);
        if length(res)==13;
          dd=[dd res(1:2:13)];
        end;
      end;
      dd=dd';
      [len(i),ncol]=size(dd);

  % save this data in a cast file
      eval([id,zer,num2str(ic),'=dd;']);
      fclose(idf);
    end

  elseif wctd==6;
%   length is not relevant here since the depths are uneven

  elseif wctd==7;
    id='el04';
    ty='.cnv';

    disp('Hardcoded to expect p,t,c,d1,d2,d3,d4,d5,s,tr in that order')
    len=zeros(ncasts,1);
    i=0;
    for jj=1:length(snum);
      if snum(jj)<10;sta='00';elseif snum(jj)>=10&snum(jj)<100;sta='0';else;sta='';end;
      fid=fopen(['c:\gbs99\ctd\el9904\',id,sta,num2str(snum(jj)),'.cnv'],'r');
      str1='*END*';
      mark=fgetl(fid);
      [mark1,rem]=strtok(mark);
      while strcmp(str1,mark1)~=1
         mark=fgetl(fid);
         [mark1,rem]=strtok(mark);
      end
      mark=fgetl(fid);
      mark=fgetl(fid);
      outfid=fopen('txtfile','w');
      while mark~=-1;
        if length(mark)>0;
           fprintf(outfid,'%s\n',mark);
        end;
        mark=fgetl(fid);        
        mark=fgetl(fid);
      end
      fclose(outfid);
      fclose(fid);
      load txtfile -ascii;
      dd=txtfile;
      %     save this data as a cast variable
      eval([id,sta,num2str(snum(jj)),'=dd;']);
      ic=snum;
      [len(jj),ncol]=size(dd);
    end
    
  end;

%%%%%%%%%%%%%%%%%%%%%%  % initialize matrices %%%%%%%%%%%%%%
  
  m=max(len);
  pmatrix=NaN*ones(m,ncasts);
  tmatrix=pmatrix;
  smatrix=pmatrix;
  dmatrix=pmatrix;
  %additional values in case of General Oceanics machines
  if wctd==1;
    rmatrix=pmatrix;
    trmatrix=pmatrix;
    fmatrix=pmatrix;
  end;
  dmatrix=pmatrix;
  sigmamatrix=pmatrix;


  %%%%%%%%%%%%%%%%%%%%%%%%%  GENERATE VARIABLE MATRICES %%%%%%%%%%%%%%%

  % fill matrices, determine min & max per cast, and save index of deepest
  % record in a variable maxpindex(i)
  for i=1:ncasts
    ic=snum(i);
    if ic<10,numz='00';elseif ic>=10&ic<100;numz='0';else;numz='';end;
    if wctd<4|wctd==5;
      if wctd==5;sing=-1;else;sing=-1;end;
      eval(['pmatrix(1:10:len(i),i)=',id,numz,num2str(ic),'(1:10:len(i),1)*sing;']);
      eval(['tmatrix(1:10:len(i),i)=',id,numz,num2str(ic),'(1:10:len(i),2);']);
      eval(['smatrix(1:10:len(i),i)=',id,numz,num2str(ic),'(1:10:len(i),3);']);
      if wctd==1;
        eval(['rmatrix(1:len(i),i)=',id,numz,num2str(ic),'(:,4);']);
      end;
    elseif wctd==4;
      idd=find(posn(:,11)==ic&posn(:,3)<=dmax*(-1));
      pmatrix(1:len(i),i)=posn(idd,3)*(-1);
      tmatrix(1:len(i),i)=posn(idd,5);
      smatrix(1:len(i),i)=posn(idd,6);
    elseif wctd==6;
%      idd=find(posn(:,2)==cumdist(i)*1000);
      test98=posn(:,2);
%      test99=-1*(cumdist(i)*1000-max(posn(:,2)));
      test99=cumdist(i)*1000;
      idd=find(test98==test99);
      pmatrix(1:len(i),i)=posn(idd,3);
%     if we want normal and tangential velocity specify col 7 & 8
%      disp('Hardcoded to pick up un and ut from rvt file') 
      tmatrix(1:len(i),i)=posn(idd,5);
      smatrix(1:len(i),i)=posn(idd,6);
%       sigmamatrix(1:len(i),i)=posn(idd,4);

    elseif wctd==7;
      eval(['pmatrix(1:len(i),i)=',id,numz,num2str(ic),'(:,3)*(-1);']);
      eval(['tmatrix(1:len(i),i)=',id,numz,num2str(ic),'(:,2);']);
      eval(['smatrix(1:len(i),i)=',id,numz,num2str(ic),'(:,3);']);
      eval(['trmatrix(1:len(i),i)=',id,numz,num2str(ic),'(:,11);']);
      eval(['fmatrix(1:len(i),i)=',id,numz,num2str(ic),'(:,7);']);
    elseif wctd==11;
      eval(['pmatrix(1:len(i),i)=-press(id(i):id(i)+len(i)-1);']);
      eval(['tmatrix(1:len(i),i)=temp_ctd(id(i):id(i)+len(i)-1);']);
      eval(['smatrix(1:len(i),i)=sal_ctd(id(i):id(i)+len(i)-1);']);
    end;
    

    sigmamatrix(1:len(i),i)=sigmat(smatrix(1:len(i),i),tmatrix(1:len(i),i))-1000; 
    mint(i)=min(tmatrix(1:len(i),i));
    maxt(i)=max(tmatrix(1:len(i),i));
    mins(i)=min(smatrix(1:len(i),i));
    maxs(i)=max(smatrix(1:len(i),i));
    if wctd==1;
      minr(i)=min(rmatrix(1:len(i),i));
      maxr(i)=max(rmatrix(1:len(i),i));
    end;
    if wctd==1|wctd==7;
      mintr(i)=min(trmatrix(1:len(i),i));
      maxtr(i)=max(trmatrix(1:len(i),i));
      minf(i)=min(fmatrix(1:len(i),i));
      maxf(i)=max(fmatrix(1:len(i),i));
    end;
    minsig(i)=min(sigmamatrix(1:len(i),i));
    maxsig(i)=max(sigmamatrix(1:len(i),i));
    [maxp(i),maxpindex(i)]=max(abs(pmatrix(1:len(i),i)));
  end

  % overall min & max
  mintemp=min(mint);
  maxtemp=max(maxt);
  minsalt=min(mins);
  maxsalt=max(maxs);
  if wctd==1;
    minrad=min(minr);
    maxrad=max(maxr);
  end;
  if wctd==1|wctd==7;
    mintran=min(mintr);
    maxtran=max(maxtr);
    minflur=min(minf);
    maxflur=max(maxf);
  end;
  minsigma=min(minsig);
  maxsigma=max(maxsig);



  % determine the max index of maximums
  [maxmax,indexmax]=max(maxpindex(1:ncasts));
  cumdist2=cumdist;

%%%%%%%%%%%%%%%%%%%%%   BLANKING POLYGON CREATED  %%%%%%%%%%%%%%%%%%%%%

%
%  write a blanking file requires getting the bottom depth somehow
%
  if vbdop==1;
  bdop=menu('How do you want to specify bottom depth?',...
          'get from posn.dat file',...
          'just add one meters to max ctd sample depth',...
          'prompt user for each station',...
          'get from a bd.dat file',...
          'read directly from a file',...
          'use beardsleys bot depth finder');
  end;

  if vdmax==1;
  % define maximum depth and xdist desired
    dmax=input('Enter Maximum Depth: ');
    if dmax>0;dmax=-dmax;end;
    vdmax=0;
  end;
  if vxmax==1;
    xmax=input('Maximum distance (km) alongtack? ');
    xmin=input('Minimum distance (km) alongtack? ');
    vxmax=0;
  end;
  disp(['Hardcoded to limit depth to ',num2str(dmax)]);
  disp(['Hardcoded to limit dist offbank to ',...
        num2str(xmin),' to ',num2str(xmax)]);

  if bdop==1; 
    bd=posn(id,3);
    a=[min(cumdist2) max(cumdist2) min(bd) dmax];
    blk=[[cumdist2(1:ncasts);cumdist2(ncasts);cumdist2(1);cumdist2(1)],...
      [bd;a(4);a(4);bd(1)]];
  elseif bdop==2;
%    a=[min(cumdist2) max(cumdist2) min(maxp) max(maxp)];
%    blk=[[xmin; (max(cumdist2)-cumdist2(1:ncasts)) ;xmax;xmax;xmin;xmin],...
%      [-maxp(1); -maxp(ncasts:-1:1)'-1;dmax;dmax;dmax;-maxp(1)]];
    blk=[[xmin; cumdist2(1:ncasts)' ;xmax;xmax;xmin;xmin],...
      [-maxp(1); -maxp(1:ncasts)'-1;dmax;dmax;dmax;-maxp(1)]];
  elseif bdop==3;
    disp(' ');
    disp('Hardcoded to get bottom dep from user');
    for i=1:ncasts;
      bd(i)=input(['Enter bottom depth for station ',num2str(snum(i)),': ']);
    end;
    a=[min(cumdist2) max(cumdist2) min(bd) max(bd)];
    blk=[[cumdist2(1:ncasts)' ;cumdist2(ncasts);cumdist2(1);cumdist2(1)],...
      [bd';a(4);a(4);bd(1)]];
  elseif bdop==4;
    disp(' ');
    disp('Hardcoded to get bottom dep from bd.dat file');
    eval(['load ',pth,'/bd.dat']);
    for i=1:ncasts;
      for kk=1:length(bd);
        if bd(kk,1)==snum(i);
          bdn(i)=bd(kk,2);
        end;
      end;
    end;
    bdn=-1*bdn;
    bd=bdn;
       
    a=[min(cumdist2) max(cumdist2) min(bdn) dmax];
    blk=[[cumdist2(1:ncasts)' ;cumdist2(ncasts);cumdist2(1);cumdist2(1)],...
      [bd';a(4);a(4);bd(1)]];
  elseif bdop==5;
    load blk.dat
  elseif bdop==6;
    
%    bd=-1*w_depth(-1*q(:,2),q(:,1));
    bd=-1*w_depth(q(:,2),q(:,1));
    bdn=bd;
    a=[min(cumdist2) max(cumdist2) min(bdn) dmax];
    blk=[[cumdist2(1:ncasts)' ;cumdist2(ncasts);cumdist2(1);cumdist2(1)],...
        [bd;a(4);a(4);bd(1)]];
  end;

% adcp skip ends here
  end;

  %%%%%%%%%%%%%%%%   NOW WE START THE CONTOURING  %%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%   DEFINE VARIABLE   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if vfield==1;
    field=menu('Please choose the field ',...
           'Temperature',...
           'Salinity',...
           'Sigma-t',...
           'transmissoelmmeter',...
           'fluorometer',...
           'Sigma-t',...
           'u',...
           'v');
  end;

% Note: we assign some default contour features while we are in this loop
  if field==1
    cfield='temp';
    zz=tmatrix;
    dincr=[8:2:22];
  elseif field==2;
    disp('Hardcoded to plot salinity')
    cfield='salt';
    zz=smatrix;
    dincr=[32:.5:36];
  elseif field==3;
    cfield='sigm';
    zz=sigmamatrix;
    dincr=[23.:.5:27];
  elseif field==4;
    cfield='trans';
    zz=trmatrix;
    dincr=[0:.2:5];
  elseif field==5;
    cfield='flur';
    zz=fmatrix;
    dincr=[0:.2:5];
  elseif field==6;
    cfield='sig';
    zz=sigmamatrix;
  elseif field>=7;
%   adcp is a special case
%   all the gridding & contouring is done in "contadcp.m"
    contadcp

%   code to show adcp from now on
    fieldn=7;

  end;
  if field<7;fieldn=1;end;

% skip most of the following if this is an adcp section
  if fieldn~=7;

% %%%%%%%%%%%%%%%%%%  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(zz);
    if incr==1;
      v=[mi:rg/5:ma];
    elseif incr==2;
      v=input('Enter min:inc:max contour int: ');
    elseif incr==3;
      disp('Hardcoded contour intervals to user defined to:');
        v=dincr;
%        v=[-10:1:10]
%       v=[4:.1:4.8 5 6 7 8]
%       v=[32.:.5:36]
    elseif incr==4;
%     attempt to intgerize the intervals
      if rg>2;
        v=floor(minnnan(zz)):floor(rg/5):ceil(maxnnan(zz));
%        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;

%   contourlabel fonsize dependent on number of sections
    if mm~=1;
      clfs=15-nsec;
    else;
      clfs=14;
    end;
    
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   GRID DATA  &&&&&&&&&&&&&&&&&&&&&&&

%    gd=menu('Grid data?','use griddata',...
%                         'use sagas mincurvi',...
%                         'use sagas kriging',...
%                         'use sagas objmap',...
%                          'none');       
    gd=4;
    if fieldn==7&ksec==1;gd=5;end;
    if gd<=4;
      ngp=25;
      disp(['Hardconlond to grid data with ',num2str(ngp),' grid pts']);
      xi=[min(cumdist):(max(cumdist)-min(cumdist))/ngp:max(cumdist)];
      yi=[minnnan(minnnan(pmatrix)):(maxnnan(maxnnan(pmatrix))-minnnan(minnnan(pmatrix)))/ngp:maxnnan(maxnnan(pmatrix))]';
      [nr,nc]=size(zz);
      x=ones(nr,nc);
      for kk=1:nc;
        x(:,kk)=cumdist(kk)*x(:,kk);
      end;
      y=pmatrix;
      
      disp('Gridding in progress ...');

%     find the Nan's and get rid of them
      xn=reshape(x,nr*nc,1);
      yn=reshape(y,nr*nc,1);
      zn=reshape(zz,nr*nc,1);
      id=find(isnan(yn)~=1);
      
%     get a mesh
      [X,Y]=meshgrid(xi,yi);
      
      if gd==1;
        zi=griddata(xn(id),yn(id),zn(id),X,Y);
      elseif gd==2;
        zi=mincurvi(xn(id),yn(id),zn(id),X,Y);
      elseif gd==3;
        zto=input('Enter some number betw 0 and 1: ');
        xysc=input('Enter vector w/x&y scales: ');
        zi=kriging(xn(id),yn(id),zn(id),X,Y,xysc,zto,[20 16 200 1/3 1]);
      elseif gd==4;
%        zto=input('Enter some number betw 0 and 1: ');
         zto=.5;
%        xysc=input('Enter vector w/x&y scales: ');
         xysc=[3 5];
        zi=objmap(xn(id),yn(id),zn(id),X,Y,xysc,zto,[40 32 200 1/3 1]);
      end;
    end;


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

    if vptype==1;
      ptype=menu('Plot type?','pcolor & contour','econtour fill','contour lines');
    end;
%%%%%%%%%%%%%%%%%%%%%%%%%%  Do the Contouring Now %%%%%%%%%%%%%%%%


      if ptype==1;
  	pcolor(cumdist2(1:ncasts) ,pmatrix(1:maxmax,innlonxmax),zz(1:maxmax,1:ncasts));
  	hold on
  	cs=contour(cumdist2(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 fill.'); 
        if gd<=4;
           cd c:\gbs99\atsee\newextc
           extcontour(xi,yi,zi,v,'fill','fontsize',clfs);
           cd c:\gbs99\atsee
        else;
  	        econtour(cumdist2(1:ncasts),pmatrix(1:maxmax,innlonxmax),zz(1:maxmax,1:ncasts),v,'fill');
        end;
           caxis([min(v) max(v)]);
         if ksec==nsec|mm==1|vfield==1;
           if vfield==1&nsec~=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;
      elseif ptype==3;

 %        path('e:\atsee\newextc',path)
        disp('Hardconlond to use EXTCONTOUR contour lines.'); 
        if gd<=4;
          cd c:\gbs99\atsee\newextc
          extcontour(xi,yi,zi,v,'label','fontsize',clfs,'color','g','linewidth',2);
          cd c:\gbs99\atsee
        else;
  	        extcontour(cumdist2(1:ncasts),pmatrix(1:maxmax,innlonxmax),zz(1:maxmax,1:ncasts),v,'fill');
        end;
      end;
 
%     heres we we might want to save this gridded field
%      for a subsequent overlay
      if vsgf==1;
        sgf=menu('Save this gridded field for subsequent overlay?','y','n');
      end;
      if sgf==1;
        xxi=xi;yyi=yi;zzi=zi;
%       save min & max yearday
%        disp(['min yrday0_gmt: ',num2str(min(yd),6),' max yrday0_gmt: ',num2str(max(yd),6)]);
%       save name according to section
        wsec=input('What section is this?');
        eval(['save sigmam',num2str(wsec),' xxi yyi zzi']);
      end;

      if wctd~=6&wctd~=9&wctd~=10;
        if ksec==1;
          for i=1:ncasts
            h=text(cumdist2(i)',0,num2str(snum(i)));
            set(h,'FontSize',8,'HorizontalAlignment','center',...
              'verticalalignment','bottom')      
          end
          if sgf==1;
%           save this info for plotting on adcp
            eval(['save sigmas',num2str(wsec),' cumdist2 ncasts snum indexmax maxp']);
          end;
        end;
      end; 
%  adcp section case can start up again
  end;


%     blank out bathymetry
      hold on;
      %      if ksec==1;
      fill(blk(:,1),blk(:,2),'k');
%      fill(blk(:,1),blk(:,2),[1 .62 .4]);
%      end;

      hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%   Fix Axis up  %%%%%%%%%%%%%%%

    if ksec==nsec|mm==1;
      if ref==1;
        xlabel('Along Track Distance (km) Relative to 1st Station')
      elseif ref==2;
        xlabel('Along Track Distance (km) Relative to Last Station')
      elseif ref==3;
        xlabel(['Distance from ',num2str(rla),'N ~ along ',num2str(rlo),'W (km)'])
      elseif ref==4;
        xlabel(['Distance Offbank 60m Isobath(km)']) 
      elseif ref==5;
        xlabel(['Distance Offbank 60m Isobath(km)'])
      end;
      ylabel('Depth (m)')
    else;
      set(gca,'xticklabels',[' ']);
    end;

%   position subplot
    if mm==2;
      set(gca,'position',[.15 .975-(ksec/nsec*0.85) .65 (1/nsec*0.75)])
    else;
      set(gca,'position',[.15 .2 .7 .6])
%      here's where we change the width of the panel depending on shelf
%      width as calculated inside "postdrft"
%      set(gca,'position',[.15*sw1/shlfwdth .2 .7*shlfwdth/sw1 .6]);
    end;

%   set axis!!
    axis([xmin xmax dmax 0]);
 
%   set colormap
    colormap(jet)


%%%%%%%%%%%%%%%%%%%%%%%%   Special text & postings   %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%        Added Lat & Lon pos'n        %%%%%%%%%%%%%%%%%%%%%
  if wctd==10;
    q=[yadcp(1,2) yadcp(1,1);yadcp(length(yadcp),2) yadcp(length(yadcp),1)];
  end;
  disp('Lat/lon labels hardcoded out');
%  text(cumdist3(1),-3,num2str(q(1,1),5))
%  text(cumdist3(length(cumdist3)),-3,num2str(q(length(q),1),5),...
%    'horizontalalignment','right')
%  text(cumdist3(1),-7,num2str(q(1,2),5))
%  text(cumdist3(length(cumdist3)),-7,num2str(q(length(q),2),5),...
%    'horizontalalignment','right')


%%%%%%%%%%%%%%%%%%%%%%%%   Special text & postings   %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%         Based on the section#    %%%%%%%%%%%%%%%%%%%%%

%     Number the section
%      text(2,-65,'Deploy#4','fontsize',14,'color','k')
%      text(2,-72,['Transect # ',num2str(ksec)],'fontsize',14,'color','k'); 

%     label the variable
%       text(0,-55,cfield,'fontsize',16,'color','k');

      if ksec==1;
       text(0,20,'CTD Cast #s');
%       text(25,-5,'Climatological Spring Sigma-t','color','k','fontsize',16);
         text(1,dmax+20,'Temperature','color','w','fontsize',16);
         title('EL9904 Transect#5 19 April 1999','fontsize',16)
         %        title(['ADCP residual ',cfield,' velocity (cm/s) & Sigma-t'],'fontsize',14)
%        title(['ADCP along-bank residual (cm/s) & Sigma-t'],'fontsize',16)
%         title(['ADCP observed Velocity (cm/s) (minus Moody @K tide +1hr)'],'fontsize',14)
%        title(['June 1997 ',cfield,' velocity (cm/s)']); 
%        text(1,-65,' w/DEPTH-AVERAGED multi-constituent Candela tide removed','fontsize',14)
%        title(['Cross-Bank ',cfield,'  Structure ']);
%        title(['Cross-Great South Channel Structure OC311 August 1997 ']);
%        title('Mooring Line CTD Temperature Structure 21 May 1997');
%        title(['Cross-Bank ',cfield,'  Structure']);
%         title('16 day particle tracks in May 1997 CTD derived 3-d velocity field',...
%          title('May 1997 ADCP detided velocity, CTD sigma-t, and  MOCNESS hauls',...
%           'fontsize',14) 
%          title('May 1997 OAX CTD-derived velocity and sigma-t w/ MOCNESS hauls',...
%           'fontsize',14) 
%         disp('click position for cbar label')
%         text(-18,-75,'along-bank (north-northeastward) velocity (cm/s) contoured',...
%           'color','k','fontsize',14)
         
      end;

     if ksec==2;
%       text(1,-65,'Sigma-t Contour Lines','fontsize',14)
        text(1,dmax+20,'Salinity','color','w','fontsize',16);
%       text(25,-5,'Climatological Spring E/W Flow','color','k','fontsize',16);
%        title(['CTD derived residual (cm/s)'],'fontsize',14)
     end;

     if ksec==3;
        text(1,dmax+20,'Fluoresence','color','w','fontsize',16);
     end;

%%%%%%%%%%%%%%%%%   Special text & postings          %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%         Based on the grid method   %%%%%%%%%%%%%%%%%%%%%


%     here we post text on gridding methods if desired
      ids=0;
      if ids==1;
      if gd==1;
        text(5,-90,'Using "griddata"','fontsize',16);
      elseif gd==2;
        text(5,-90,'Using "mincurvi"','fontsize',16);
      elseif gd==3;
        text(5,-90,'Using "kriging"','fontsize',16);
      elseif gd==4;
        text(5,-90,'Using "objmap"','fontsize',16);
      elseif gd==5;
        text(5,-90,'Using no gridding','fontsize',16);
      end;
      end;

%%%%%%%%%%%%%%%%%   Special text & postings          %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%         Based on data type (wctd)  %%%%%%%%%%%%%%%%%%%%%

%   add date text (get date depending on format of input
    if wctd==4;
%      disp(['Hardcoded date to ',yr]);
%      dd=calday(str2num(yr),posn(idd(1),4));
%      ddm=caldaym(dd(2));
%      text(3,dmax+7,[num2str(dd(3)),num2str(ddm),' 19',yr],'fontsize',14);
%      text(30,-10,[num2str(dd(3)),num2str(ddm),' 19',yr],...
%         'color','k','fontsize',14);
%        yadcp=[q(1,2) q(1,1)];
%        postmoc
%        miny=131;maxy=148;yadcp=[-67.7795 41.1469];
%        postdrft
        pmoor=1;
        if pmoor==1;
          postmoor;
        end;
    elseif wctd==5;
      disp(['Hardcoded date to 19',num2str(yr)]);
      ddd=calday(yr,posn(id2(1),2));
      dddm=caldaym(ddd(2));
      text(30,-10,[num2str(ddd(3)),num2str(dddm),' 19',num2str(yr)],...
        'color','k','fontsize',14);
    elseif wctd==6;
%       title([num2str(fd),' May 1997 at ',num2str(fh),'hrs (UTC)'],'fontsize',16);

%     Here is where we post velocity arrows on the section
%     with a combination of "sizesym" and then "vecplot2"

       hold on;
%       load /data5/jmanning/model/g2s.5b/vel/MJ/resid4.nrv
%       vel=resid4;
       vel=posn;
%       subsample vel 
        ssr=1;

% note: for .rvt files use col 8 & 9 for v & w
        pause(1);
        disp('vecplot');
        vecplot2(vel(1:ssr:length(vel),2)/1000,vel(1:ssr:length(vel),3),...
         vel((1:ssr:length(vel)),6),1000*vel((1:ssr:length(vel)),7),...
         5,'cm/s',-3,-70,'r')
        disp('vecplot');

%        now post drifters on this
%        routine assume miny and maxy is a time range of interest
%        and yadcp matrix has long & lat in col 1 & 2
%        miny=ydnow-(2/24);maxy=ydnow;yadcp=[-67.7795 41.1469];
%        disp('overlays <.25 hours of drifters');
         miny=131;maxy=148;yadcp=[-67.7795 41.1469];
%         postdrft;

        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.22:.1:25.82],'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;
%        yadcp=[-67.850   40.950];       
%        postmoc
         
    elseif wctd==9;

%     posting done within "contadcp"
         
    elseif wctd==11;
      ddd=calday(yr,meannnan(yrday0_gmt));
      dddm=caldaym(ddd(2));
      text(30,-10,[num2str(ddd(3)),num2str(dddm),' 19',num2str(yr)],...
          'color','k','fontsize',14);
    end;


%%%%%%%%%%%%%%%%%%%%%%%%   CAPTURE THIS IMAGE

    if mm==1;
        disp(['Saving gif file... ']);
        eval(['print -dgif8 v',num2str(fd),num2str(fh-1)]);
%       add to list file
        fprintf(fid2,'%s\n',['v',num2str(fd),num2str(fh-1),'.gif']);
    end;
    disp('Any Annotations?');
    pause(1);

%   clear q  blk snum cumdist2 cumdist bd bdn;
    clear cumdist
  end
end
    orient tall;
if mm==1;
  disp('go to whsun1 and type:')
  disp(' ~/flic/xanim2703/fbm2fli list.dat this.flc');
  disp(' and then on quoody type: ')
  disp('xanim this.flc');
end;
    


