%Thermohaline Plots
% WRITTEN by Theo Touvet in July 2010
clf
colormap(mk_colmap)

%Plot Temperature
PT = Tsw;
PT(1:Nbcs)=NaN;
subplot(4,1,1)
surf(XP,YP,PT(NodeP(2:app.Ny+1,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('Tsw(^oC)');
TT(4)=title(sprintf('Temperature (^oC), (app.Nx=%g, app.Ny=%g)',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp

%Plot Salinity
PS = Ssw;
PS(1:Nbcs)=NaN;
subplot(4,1,2)
surf(XP,YP,PS(NodeP(2:app.Ny+1,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('S(ppm)');
TT(4)=title(sprintf('Salinity (psu)'));
formatplot([],gca,TT);view(2),colorbar,shading interp

%Plot Density
% NON-LINEAR equation of state using SW_DENS0: denisty of sea water AT
% ATMOSPHERIC PRESSURE
rho_0=1028; g=9.81;
rho = g*(sw_dens0(Ssw,Tsw)-rho_0*ones(size(P)))/rho_0; %perturbation density divided by rho_0 and multiplied by g. Here S has to be in psu and Tsw in Celsius degree!
Prho = rho;
Prho(1:Nbcs)=NaN;
subplot(4,1,3)
surf(XP,YP,Prho(NodeP(2:app.Ny+1,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('rho(kg/m^3)');
TT(4)=title(sprintf('"Density" (m^2/s)'));
formatplot([],gca,TT);view(2),colorbar,shading interp

%Quiver Plot
subplot(4,1,4)
tmp1 = 0.5*(u(Nodeu(2:app.Ny+1,1:app.Nx))+u(Nodeu(2:app.Ny+1,2:app.Nx+1)));
tmp2 = 0.5*(v(Nodev(2:app.Ny+1,2:app.Nx+1))+v(Nodev(1:app.Ny,2:app.Nx+1)));

stepsize = 3; % useful to change the number of arrows in the "quiver" plot (the smaller "stepsize" is, the more you have arrows)
quiver(XP(1:stepsize:size(XP,1),1:stepsize:size(XP,2)),YP(1:stepsize:size(YP,1),1:stepsize:size(YP,2)),(tmp1(1:stepsize:size(tmp1,1),1:stepsize:size(tmp1,2))),(tmp2(1:stepsize:size(tmp2,1),1:stepsize:size(tmp2,2))),2,'g') 
hold on
streamline(XP,YP,tmp1,tmp2,[-a/2+dx, a/2-dx, -a/4, a/4],[b/2, b/2, b/2, b/2])

surf(XP,YP,sqrt(tmp1.^2+tmp2.^2),'FaceAlpha',0.2,'EdgeAlpha',0.0), colorbar
hold off
xlim([-a/2 a/2]), ylim([0 b])
TT(1)=xlabel('x(m)');
TT(2)=ylabel('y(m)');
TT(3)=zlabel('V(m/s)');
TT(4)=title(sprintf('Velocity, time =%03.03g',T));

formatplot([],gca,TT);view(2)
