%% EPS131, Eli Tziperman. 2016. 
%% Demonstrating: 
%% (1) Cabeling
%% (2) Sigma theta is not appropriate for stability,
addpath('~/matlab/seawater_eqn_of_state_gsw_matlab_v3_06/' ...
        ,'~/matlab/cbfreezColors-cm_and_cb_utilities/')
set(0,'defaulttextfontsize',14); set(0,'defaultaxesfontsize',14);

%% Demonstrate cabeling:
%%%%%%%%%%%%%%%%%%%%%%%%
S1=34:0.01:36;
T1=-2:0.05:30;
[T,S]=meshgrid(T1,S1);
p0=0;
p4=4000;
sigma0=gsw_rho(S,T,p0)-1000;
sigma4=gsw_rho(S,T,p4)-1000;
%  S  =  Absolute Salinity                      [ g/kg ]
%  T  =  Conservative (~ potential) Temperature [ deg C ]
%  p  =  pressure                               [ dbar ]

%% salinity and temperature to demonstrate cabeling:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Two points demonstrating cabeling:
S1=34.99441;
T1=-0.6;
S2=35.64666;
T2=5.45;
% Two points demonstrating exchange of stability between sigma_0 and sigma_4:
S3=35.61;
T3=4.4;
S4=34.92;
T4=-0.2;
T_mean=(T1+T2)/2;
S_mean=(S1+S2)/2;


%% contour plots of density at p=0 and 4000:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1);clf
subplot(1,2,1);
contourf(S,T,sigma0,'edgecolor','none')
jet_inv=jet; jet_inv=jet_inv(end:-1:1,:); colormap(jet_inv);
h=colorbar;
caxis([min(sigma0(:)),max(sigma0(:))]);
hold on
[C,hc]=contour(S,T,sigma4,'k'); clabel(C,hc,'fontsize',14,'LabelSpacing',400);
cbfreeze(h);
title('color: \sigma_0,   black: \sigma_4')
xlabel('S')
ylabel('T')

%% add to the contour plot the temp and sal of the two points whose
%% mixing leading to cabeling:
hold on
plot(S1,T1,'yx','markersize',14,'linewidth',2)
plot(S2,T2,'yx','markersize',14,'linewidth',2)
plot(S_mean,T_mean,'r+','markersize',14,'linewidth',2)
line([S1,S2],[T1,T2]);
text(34.1,3,'cabeling','color','w' ...
     ,'fontsize',16)

subplot(1,2,2);
contourf(S,T,sigma4,'edgecolor','none')
jet_inv=jet; jet_inv=jet_inv(end:-1:1,:);colormap(jet_inv);
h=colorbar;
caxis([min(sigma4(:)),max(sigma4(:))]);
hold on
[C,hc]=contour(S,T,sigma0,'k'); clabel(C,hc,'fontsize',14,'LabelSpacing',400);
cbfreeze(h)
title('color: \sigma_4,   black: \sigma_0')
xlabel('S')
ylabel('T')

%% add to the contour plot the temp and sal of the two points whose
%% stability changes between sigma_0 and sigma_4000:
hold on
plot(S3,T3,'wx','markersize',14,'linewidth',2)
plot(S4,T4,'wx','markersize',14,'linewidth',2)
text(34.1,3,'stability: \sigma_0 vs \sigma_4','color','w' ...
     ,'fontsize',16)

fprintf(1,['\n*** Mixing two water masses with the same\n' ...
           '    density leads to denser water (cabeling):\n']);
sigma1=gsw_rho(S1,T1,0)-1000;sigma2=gsw_rho(S2,T2,0)-1000;
sigma_mean=gsw_rho(S_mean,T_mean,0)-1000;
fprintf(1,'T1=%g, S1=%g, sigma(S1,T1,0)=%7.5f;\n',T1,S1,sigma1);
fprintf(1,'T2=%g, S2=%g, sigma(S2,T2,0)=%7.5f;\n',T2,S2,sigma2);
fprintf(1,'T_mean=%g, S_mean=%g, sigma(S_mean,T_mean,0)=%7.5f.\n' ...
        ,T_mean,S_mean,sigma_mean);

%% Stability and sigma theta:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S3=34.99441;
T3=-0.6;
p1=0;
S4=35.8;
T4=5.45;
p2=4000;
sigma1p1=gsw_rho(S3,T3,p1)-1000;
sigma2p1=gsw_rho(S4,T4,p1)-1000;
sigma1p2=gsw_rho(S3,T3,p2)-1000;
sigma2p2=gsw_rho(S4,T4,p2)-1000;
fprintf(1,['\n*** Water mass S3,T3 is lighter than S4,T4\n' ...
           '      at the surface but denser at depth:\n']);
fprintf(1,'T3=%g, S3=%g, p1=%g; T4=%g, S4=%g, p2=%g;\n' ...
        ,T3,S3,p1,T4,S4,p2);
fprintf(1,' sigma(S3,T3,p1)=%7.5f < sigma(S4,T4,p1)=%7.5f,\n' ...
        ,sigma1p1,sigma2p1);
fprintf(1,' sigma(S3,T3,p2)=%7.5f > sigma(S4,T4,p2)=%7.5f,\n' ...
        ,sigma1p2,sigma2p2);

if 0
  %% save as pdf:
  set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [10 5]);
  set(gcf, 'PaperPosition', [0 0 10 5]); % [left, bottom, width, height];
  saveas(gcf,sprintf('Figures/Figure-nonlinear-eqn-of-state.pdf'));
end
