% plot T, to understand it
% M is the number of nodes,
% N is the number of triangles

clf
clear all
load mesh

filehist='gom_history.cdf';
filedir='/home/pringle/fvcom/fvcom11jul'
[jnk,basedir]=system(['basename ' filedir]); basedir=basedir(1:end-1)
nchist=netcdf([filedir '/outmedm/' filehist],'nowrite');
t=nchist{'t1'};
s=nchist{'s1'};
u=nchist{'u'};
v=nchist{'v'};
thour=nchist{'thour'};


%define the time and depth to get
 zplot=-20
 %nt=237;
 %nt=length(thour)-8
 nt=5*4-1

dnt=8; %data written out every 8 tidal cycles. 

%get that chunk
tchunk=squeeze(t(nt:(nt+dnt-1),:,:));
schunk=squeeze(s(nt:(nt+dnt-1),:,:));
uchunk=squeeze(u(nt:(nt+dnt-1),:,:));
vchunk=squeeze(v(nt:(nt+dnt-1),:,:));
time=thour(nt:(nt+dnt-1));

%take time mean from nt:(nt+dnt-1)
time=squeeze(mean(time));
tchunk=squeeze(mean(tchunk,1));
schunk=squeeze(mean(schunk,1));
uchunk=squeeze(mean(uchunk,1));
vchunk=squeeze(mean(vchunk,1));

%interp t to depth
 tic;
 tdepth=nan*ones(size(squeeze(tchunk(1,:))));
 sdepth=nan*ones(size(squeeze(schunk(1,:))));
 for m=1:M
   %if (rem(m,1000)==0)
   %  m
   %end
   tdepth(m)=interp1_fast(mesh.nodez(:,m),tchunk(:,m),zplot);
   sdepth(m)=interp1_fast(mesh.nodez(:,m),schunk(:,m),zplot);
 end
 toc

 %interp u and v to depth, but only on points to be ploted. 
 load minspace10000
 gp=goodpts;
 ud=zeros(length(gp),1);
 vd=zeros(length(gp),1);
 tic
 for g=1:length(gp)
   %if (rem(g,1000)==0)
   %  g
   %end
   ud(g)=interp1_fast(mesh.zuv(:,gp(g)),uchunk(:,gp(g)),zplot);
   vd(g)=interp1_fast(mesh.zuv(:,gp(g)),vchunk(:,gp(g)),zplot);
 end
 toc
 
 
 %draw coastlines
 col=0.7*[1 1 1];
 jnk=patch('Vertices',mesh.nodexy/1e3,'Faces',mesh.trinodes,...
     'FaceColor',col,'EdgeColor',col);

%compute density -- note, note useing Chen's eq of state, but unesco
%THIS IS A FLAW
 pres=sw_pres(-zplot,42);
 dens_depth=sw_dens(sdepth,tdepth,pres);
 sigma_depth=dens_depth-1000;
 
%plot vertices
 hold on
 patch('Vertices',mesh.nodexy/1e3,'Faces',mesh.trinodes,'Cdata',sigma_depth,...
     'edgecolor','none','facecolor','interp')
 hold off
 
 xlabel(sprintf(...
     'tidal mean density at z=%2.0f for file %s at time %4.1f hours',...
     zplot,basedir,time))

 ax=[505 1775 -660 355]; %whole domain
 %ax=(1.0e+03 *[0.6659    1.7782   -0.4752    0.3990]);
 %ax=[850 1350 -300 250];
 %ax=[1110 1325 -144 40];
 ax=[0.7477    1.5550   -0.3348    0.3371]*1e3;
 axis(ax)
 %caxis([2 25])
 axis equal

 
 %draw arrows
 %the number of seconds so that the length of the arrow = dts*U
 dts=8.64e4*1.0;

 %load what points to draw arrows on
 %load allspace
 gp=goodpts;
 
 dx=dts*ud';
 dy=dts*vd';
 offset=[dx;dy]';

 %don't plot arrows where their are no currents, or where they are less
 %than min_plot
 min_plot=0.01;
 bp=find((~isnan(ud)).*(min_plot<sqrt(ud.^2+vd.^2)).*...
     (mesh.uvnode(gp,1)>ax(1)*1e3).*...
     (mesh.uvnode(gp,1)<ax(2)*1e3).*...
     (mesh.uvnode(gp,2)>ax(3)*1e3).*...
     (mesh.uvnode(gp,2)<ax(4)*1e3));
 
 axis(ax)
 arrow(mesh.uvnode(gp(bp),:)/1e3,(mesh.uvnode(gp(bp),:)+offset(bp,:))/1e3,...
     'length',2,'tipangle',20,'width',0.25);
 

 title(sprintf(...
    'file=%s, arrow is distance traveled in %2.1f days,min arrow speed=%4.3fm/s',...
     basedir,dts/8.64e4,min_plot))

 %plot a scale arrow.
  xa=ax(1)+0.1*(ax(2)-ax(1));
  ya=ax(3)+0.9*(ax(4)-ax(3));
  ya_txt=ax(3)+0.93*(ax(4)-ax(3));
 
  arrow([xa ya],[xa+dts*0.3/1e3 ya],...
     'length',4,'tipangle',20,'width',0.5);
  text(xa,ya_txt,'30 cm/s') 
 
  axis(ax)

 
 %draw depth contours 
 load depthgrd
 hold on
 contour(depthgrd.xvec/1e3,depthgrd.yvec/1e3,...
     depthgrd.dpth,-[-270 -200 -100:20:0],'k-')
 hold off

 set(gcf,'renderer','painters')
 colorbar


 
%clean up