%
% Load data from files
%
   clear all
   !ls *.m3d *.o3d
   filename=input('Enter the name of observational file: ','s');
   [obs,gridname,year,ncol,gmt]=read_obs(filename);
%
% Filter data
%
   size(obs);
   nnv=21;
   nn=ans(1)/21;
   ii=0;
   for i=1:10:nn
      for j=1:nnv
         ij=(i-1)*nnv+j;
         ii=ii+1;
         obs1(ii,:)=obs(ij,:);
      end
   end
   clear obs;
   obs=obs1;
%
% Enter mesh generation info
%
   for i=1:ncol
      fprintf(1,' first col %4.0f entry is %12.4e \n',i,obs(1,i))
   end
   xcol=input('Enter the column number which corresponds to the x axis: ');
   ycol=input('Enter the column number which corresponds to the y axis: ');
   x=obs(:,xcol);
   y=obs(:,ycol);
%
% Create mesh and create boundary
%
   nnv=21;
   size(x);
   nn=ans(1)
   nseg=nn/nnv-1
   femgen
   bnd=detbndy(in);
%
% Begin plotting loop
%
newplot='y';
while newplot == 'y',
%
% Plot boundary
%
   figure
   whitebg('w')
   hold on
   bndo=plotbnd(x,y,bnd);
   set(bndo,'Color','k')
   axis('off')
%
% Set contour levels and make plot
%
   newcont='y';
   while newcont == 'y',
      for i=1:ncol
         fprintf(1,' first col %4.0f entry is %12.4e \n',i,obs(1,i))
      end
      scol=input('Enter the column number which corresponds to the scalar you desire to contour: ');
      scalar=obs(:,scol);
      scrange(scalar)
      cint=input('Enter the contour interval: ');
      cmin=cint*ceil(min(scalar)/cint);
      cmax=cint*floor(max(scalar)/cint);
      clear cval
      i=1;
      cval(i)=cmin;
      while cval(i) < cmax
         i=i+1;
         cval(i)=cval(i-1)+cint;
      end
      cval
      hc=lcontour2(in,x,y,scalar,cval);
      title(filename)    
      zoom on
      newcont=input('New contour interval? (y/n): ','s');
      if newcont == 'y',
         delete(hc);
      end
   end
   newplot=input('New plot? (y/n): ','s');
end