%==========================================================================
% FILE: zenodo_tailleux_os_2021_figures
%
% DESCRIPTION:
%   Matlab script for creating and plotting figures of the paper:
%   Tailleux, R., 2021: Spiciness theory revisited, with new views on
%   neutral density, orthogonality and passiveness. Ocean Science, 17,
%   203-217. 
%   https://doi.org/10.5194/os-17-203-2021
%
% DEPENDENCIES:
%   This script calls the following subroutines:
%   - set_graphics_properties
%   - get_woce_params
%   - grad_density
%   - grad_locally_referenced
%   - sin_angle
%   - gsw_gammat_analytic_os_2021
%   - gsw_pspi
%   - plot_figure_correlations
%   - plot_figure_histograms_angle
%   - plot_figure_introduction_revised
%   - spice_tailleux
%
% EXTERNAL LIBRARIES:
%   - Gibbs Sea Water library available at: www.teos-10.org
%   - cmocean library available at: 
%   https://uk.mathworks.com/matlabcentral/fileexchange/57773-cmocean-perceptually-uniform-colormaps
%
% DATASETS:
%   - waghc_params.nc: available at:
%   https://www.cen.uni-hamburg.de/icdc/data/ocean/woce-climatology.html
%   - gamman_woce.mat 
%   Matlab file storing the value of neutral density independently computed
%   by using the neutral density software available at:
%   http://www.teos-10.org/preteos10_software/neutral_density.html
%
%
% AUTHOR:
%   Remi Tailleux, University of Reading, 12 June 2021 
%
% DISCLAIMER:
%   This software is offered as is and without warranty to be bug-free. If
%   you find any bug/problem and/or are able to propose any improvement to
%   my code, I'd more than happy to hear from you. In that case, please 
%   email me at: R.G.J.Tailleux@reading.ac.uk 
%
%   If you make significant use of the code for your research, please
%   acknowledge this by citing my Ocean Science paper above. 
%   Tailleux, R. 2021: Spiciness theory revisited, with new views on
%   neutrality, orthogonality and passiveness. Ocean Science, 17, 203-217. 
%
%==========================================================================
%% set graphics properties
%-------------------------
set_graphics_properties

%% ----------------------------------------
[longitude, latitude, zlev, ...
    bot_dep, pressure, temp, tpot, salt, ...
    gamman, sigma0, sigma2, sigma4, valid] = get_woce_params;

% Compute dimensions of three-dimensional arrays
[npi,npj,npk] = size(gamman);

% Define Earth radius
ra = 6371.d3;

disp('WOCE data read');

%% Compute three-dimensional versions of longitude, latitude, depth
%------------------------------------------------------------------
zlev3d = zeros(npi,npj,npk);
for k=1:npk
    zlev3d(:,:,k) = zlev(k);
end
lat3d = zeros(npi,npj,npk);
for j=1:npj
    lat3d(:,j,:) = latitude(j);
end
lon3d = zeros(npi,npj,npk);
for i=1:npi
    lon3d(i,:,:) = longitude(i);
end

%% Regress pressure against depth to determine conversion factor
%---------------------------------------------------------------
b = regress(pressure(valid),zlev3d(valid));

%% Concert WOCE practical salinity to Reference composition salinity
sr = salt;
sr(valid) = gsw_SR_from_SP(salt(valid)); 
sr(~valid) = NaN;

%% Compute potential temperature using TEOS10 routine 
pt0 = zeros(npi,npj,npk);
pt0(valid) = gsw_pt0_from_t(sr(valid),temp(valid),pressure(valid));
pt0(~valid) = NaN;

%% Compute COnservative Temperature
ct = zeros(npi,npj,npk);
ct(valid) = gsw_CT_from_t(sr(valid),temp(valid),pressure(valid));
ct(~valid) = NaN;

%% Computation of analytic thermodynamic neutral density gammaT_analytic
[zbis,pbis,sigbis,gammat]= gsw_gammat_analytic_os_2021(sr,ct);

%% Compute Potential densities using TEOS10 routines
sigma0(valid) = gsw_rho_CT_exact(sr(valid),ct(valid),0.)-1000.;
sigma0(~valid) = NaN;
sigma2(valid) = gsw_rho_CT_exact(sr(valid),ct(valid),2000.)-1000.;
sigma2(~valid) = NaN;
sigma4(valid) = gsw_rho_CT_exact(sr(valid),ct(valid),4000.)-1000.;
sigma4(~valid) = NaN; 
sigma1 = gsw_rho_CT_exact(sr,ct,1000.)-1000.; sigma1(~valid) = NaN; 

%% Load neutral densities file
%load('gamman_woce.mat')
gamman_jmd = gamman; % This uses the WOCE gamman. Alternatively, use the load command to use independently recomputed gamman

%% Compute neutral density vector 
gradn_rho = zeros(npi,npj,npk,3); 
[gradn_rho(:,:,:,1),gradn_rho(:,:,:,2),gradn_rho(:,:,:,3)] = ...
    grad_locally_referenced(sr,ct,pressure,longitude,latitude,zlev,ra); 
drho_dz = gradn_rho(:,:,:,3); 


%% Compute the fictitious diffusivity for gammat analytical to demonstrate 
% usefulness as compared for sigma1 
gradn_gt = zeros(npi,npj,npk);
[gradn_gt(:,:,:,1), gradn_gt(:,:,:,2), gradn_gt(:,:,:,3)] = ...
        grad_density(gammat,longitude,latitude,zlev,ra);
sin_gtbis = sin_angle(gradn_rho, gradn_gt); 
ficti_gtbis = 1000.*sin_gtbis.^2; 

gradn_gn = zeros(npi,npj,npk);
[gradn_gn(:,:,:,1), gradn_gn(:,:,:,2), gradn_gn(:,:,:,3)] = ...
        grad_density(gamman_jmd,longitude,latitude,zlev,ra);
sin_gn = sin_angle(gradn_rho, gradn_gn); 
ficti_gn = 1000.*sin_gn.^2; 

gradn_s2 = zeros(npi,npj,npk);
[gradn_s2(:,:,:,1), gradn_s2(:,:,:,2), gradn_s2(:,:,:,3)] = ...
        grad_density(sigma2,longitude,latitude,zlev,ra);
sin_s2 = sin_angle(gradn_rho, gradn_s2); 
ficti_s2 = 1000.*sin_s2.^2; 

gradn_s1 = zeros(npi,npj,npk);
[gradn_s1(:,:,:,1), gradn_s1(:,:,:,2), gradn_s1(:,:,:,3)] = ...
        grad_density(sigma1,longitude,latitude,zlev,ra);
sin_s1 = sin_angle(gradn_rho, gradn_s1); 
ficti_s1 = 1000.*sin_s1.^2; 

gradn_s4 = zeros(npi,npj,npk);
[gradn_s4(:,:,:,1), gradn_s4(:,:,:,2), gradn_s4(:,:,:,3)] = ...
        grad_density(sigma4,longitude,latitude,zlev,ra);
sin_s4 = sin_angle(gradn_rho, gradn_s4); 
ficti_s4 = 1000.*sin_s4.^2; 

gradn_s0 = zeros(npi,npj,npk);
[gradn_s0(:,:,:,1), gradn_s0(:,:,:,2), gradn_s0(:,:,:,3)] = ...
        grad_density(sigma0,longitude,latitude,zlev,ra);
sin_s0 = sin_angle(gradn_rho, gradn_s0); 
ficti_s0 = 1000.*sin_s0.^2; 

%% Figure with shaded topographic area
ficti_s0(ficti_s0>1.d-5) = NaN;
ficti_s1(ficti_s1>1.d-5) = NaN;
ficti_s2(ficti_s2>1.d-5) = NaN;
ficti_s4(ficti_s4>1.d-5) = NaN;
ficti_gtbis(ficti_gtbis>1.d-5) = NaN; 
ficti_gn(ficti_gn>1.d-5) = NaN;

%% Arrays for filling topography
ztop = min(0,-squeeze(bot_dep(661,:,1)));
ztop = max(ztop,-5500)';
zbot = -ztop*0.-6000';
ztop(1)=-6000.;ztop(end)=-6000.;
latlat = [latitude,fliplr(latitude)];
ztopzlow = [zbot,fliplr(ztop)];
latlat = [-80,-80;latlat];
ztopzlow = [-6000.,-6000.;ztopzlow];
ztopzlow(2,:) = [-6000.,0.];

%% Multipanel figure for paper 
figure
subplot(2,2,1)
pcolor(latitude,-zlev,squeeze(log10(ficti_gn(661,:,:)))');shading flat;
colorbar; cmocean('dense')
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)');
title(['(a) K_f( \gamma^n ) - 30' char(176) ' W'])
caxis([-10 -5])
axis([-80 60 -6000 0])
%pbaspect([1.1 1 1]) 
set(gca,'FontSize',10)
hold on
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

subplot(2,2,2)
pcolor(latitude,-zlev,squeeze(log10(ficti_gtbis(661,:,:)))');shading flat;
colorbar; cmocean('dense')
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)');
title(['(b) K_f( \gamma^T_{analytic} ) - 30' char(176) ' W'])
caxis([-10 -5])
axis([-80 60 -6000 0])
%pbaspect([1.1 1 1])
set(gca,'FontSize',10)
hold on
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

subplot(2,2,3)
pcolor(latitude,-zlev,squeeze(log10(ficti_s1(661,:,:)))');shading flat;
colorbar; cmocean('dense')
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)');
title(['(c) K_f( \sigma_1 ) - 30' char(176) ' W'])
caxis([-10 -5])
axis([-80 60 -6000 0])
%pbaspect([1.1 1 1])
set(gca,'FontSize',10)
hold on
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

subplot(2,2,4)
pcolor(latitude,-zlev,squeeze(log10(ficti_s2(661,:,:)))');shading flat;
colorbar; cmocean('dense')
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)');
title(['(d) K_f( \sigma_2 ) - 30' char(176) ' W'])
caxis([-10 -5])
axis([-80 60 -6000 0])
%pbaspect([1.1 1 1])
set(gca,'FontSize',10)
hold on
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

print('fig04_10pt.eps','-depsc2') 


%% Compare gammat analytical bis with gamman 
figure
pcolor(latitude,-zlev,squeeze(gamman_jmd(661,:,:))');shading flat;colorbar
caxis([25 28.5])
axis([-80 60 -6000 0])
cmocean('dense')
hold on
%contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))',27:0.1:30,'w', ...
%    'ShowText','on')
%contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))',26:0.2:26.8,'w', ...
%    'ShowText','on') 
[C2,h2]=contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))', ...
    [27.4:0.1:29],'LineColor','w');
clabel(C2,h2,'Color','w')
[C3,h3]=contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))',...
    [26:0.2:27.2],'w');
clabel(C3,h3,'Color','w')
contour(latitude,-zlev,squeeze(gammat(661,:,:))',27.4:0.1:29,'y')
contour(latitude,-zlev,squeeze(gammat(661,:,:))',26:0.2:27.2,'y') 
xlabel('Latitude (deg)');ylabel('Depth (m)')
title(['(b) \gamma^n versus \gamma^T_{analytic} - 30' char(176) ' W']) 
fill(latlat,ztopzlow,[0.8 0.8 0.8])
hold off
print('gammat_versus_gamman_bis_30W.eps','-depsc2') 

%% Reference pressure
figure
pcolor(latitude,-zlev,squeeze(pbis(661,:,:))');shading flat;colorbar
caxis([0 6000])
cmocean('deep')
axis([-80 60 -6000 0])
hold on
%contour(latitude,-zlev,squeeze(pref(661,:,:))',0:500:2000,'k', ...
%    'ShowText','on')
[C3,h3]=contour(latitude,-zlev,squeeze(pbis(661,:,:))',...
    [0:500:6000],'Color','m');
clabel(C3,h3,'Color','m')
xlabel('Latitude (deg)');ylabel('Depth (m)')
fill(latlat,ztopzlow,[0.8 0.8 0.8])
title(['(a) Reference pressure (dbar) - 30' char(176) ' W']) 
hold off
print('ref_pressure_bis_30W.eps','-depsc2') 

%% Combine plots
figure
subplot(2,1,1)
pcolor(latitude,-zlev,squeeze(pbis(661,:,:))');shading flat;colorbar
caxis([0 6000])
cmocean('deep')
axis([-80 60 -6000 0])
hold on
%contour(latitude,-zlev,squeeze(pref(661,:,:))',0:500:2000,'k', ...
%    'ShowText','on')
[C3,h3]=contour(latitude,-zlev,squeeze(pbis(661,:,:))',...
    [0:500:6000],'Color','m');
clabel(C3,h3,'Color','m','FontSize',8)
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.8 0.8 0.8])
%pbaspect([2.4 1 1])
title(['(a) Reference pressure (dbar) - 30' char(176) ' W']) 
hold off
subplot(2,1,2)
pcolor(latitude,-zlev,squeeze(gamman_jmd(661,:,:))');shading flat;colorbar
caxis([25 28.5])
axis([-80 60 -6000 0])
cmocean('dense')
hold on
%contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))',27:0.1:30,'w', ...
%    'ShowText','on')
%contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))',26:0.2:26.8,'w', ...
%    'ShowText','on') 
[C2,h2]=contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))', ...
    [27.4:0.1:29],'LineColor','w');
clabel(C2,h2,'Color','w','FontSize',8)
[C3,h3]=contour(latitude,-zlev,squeeze(gamman_jmd(661,:,:))',...
    [26:0.2:27.2],'w');
clabel(C3,h3,'Color','w','FontSize',8)
contour(latitude,-zlev,squeeze(gammat(661,:,:))',27.4:0.1:29,'y')
contour(latitude,-zlev,squeeze(gammat(661,:,:))',26:0.2:27.2,'y') 
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(b) \gamma^n versus \gamma^T_{analytic} - 30' char(176) ' W']) 
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.8 0.8 0.8])
%pbaspect([2.5 1 1])
hold off
print('fig03_12pt.eps','-depsc2')

%% Compute spiciness variables
tau_huang = gsw_pspi(sr,ct,1000.); tau_huang(~valid) = NaN;
tau_jmd = gsw_spiciness1(sr,ct); tau_jmd(~valid) = NaN;
tau_remi = spice_tailleux(sr,ct,1000.,35.,20.); tau_remi(~valid) = NaN;
tau_ql1 = spice_tailleux(sr,ct,1000.,35.,20.);
tau_ql = spice_tailleux(sr,ct,pressure,35.,20.);
%tau_remi_ref = spice_tailleux(sr,ct,pref,35.,20.); tau_remi(~valid) = NaN;
tau_remi_ref = spice_tailleux(sr,ct,pbis,35.,20.); 
tau_remi_ref(~valid) = NaN;  


%% Compute patched potential spicity
pi_ppd = zeros(npi,npj,npk);
phuang = [0.,500.,1000.,2000.,3000.,4000.,5000.];
for k=1:npk
    switch zlev(k)
        case zlev(k) < 260.
            pi_ppd(:,:,k) = gsw_pspi(sr(:,:,k),ct(:,:,k),0.);
        case zlev(k) >= 260. & zlev(k) < 760.
            pi_ppd(:,:,k) = gsw_pspi(sr(:,:,k),ct(:,:,k),500.);
        case zlev(k) >= 760. & zlev(k) < 1510.
            pi_ppd(:,:,k) = gsw_pspi(sr(:,:,k),ct(:,:,k),1000.);
        case zlev(k) >=1510. & zlev(k) < 2510.
            pi_ppd(:,:,k) = gsw_pspi(sr(:,:,k),ct(:,:,k),2000.);
        case zlev(k) >= 2510. & zlev(k) < 3510.
            pi_ppd(:,:,k) = gsw_pspi(sr(:,:,k),ct(:,:,k),3000.);
        case zlev(k) >= 3510. & zlev(k) < 4510.
            pi_ppd(:,:,k) = gsw_pspi(sr(:,:,k),ct(:,:,k),4000.);
        otherwise 
            pi_ppd(:,:,k) = gsw_pspi(sr(:,:,k),ct(:,:,k),5000.);
    end
end 
disp('computation of patched potential pi completed') 

%% Compute patched potential spiciness from McDougall stuff
tau_ppd = zeros(npi,npj,npk);
for k=1:npk
    switch zlev(k)
        case zlev(k) < 1000.
            tau_ppd(:,:,k) = gsw_spiciness0(sr(:,:,k),ct(:,:,k));
        case zlev(k) >= 1000. & zlev(k) < 2000. 
            tau_ppd(:,:,k) = gsw_spiciness1(sr(:,:,k),ct(:,:,k));
        otherwise 
            tau_ppd(:,:,k) = gsw_spiciness2(sr(:,:,k),ct(:,:,k));
    end
end
disp('computation of patched potential spiciness completed') 

%% Compute Huang's spicity variable using interpolation between pressure
% levels 

pi_ref = zeros(npi,npj,npk);
pi_ref(:) = NaN; 
phuang = [0.,500.,1000.,2000.,3000.,4000.,5000.];

validpi = (sr>10.&sr<40.) & (ct>-2.&ct<40.);

disp('Compute Huang spicity') 
for k=1:npk
    disp(k) 
    for j=1:npj
        for i=1:npi
            if validpi(i,j,k)
                for l=1:7
                    y(l) = gsw_pspi(sr(i,j,k),ct(i,j,k),phuang(l));
                end
                pi_ref(i,j,k)=interp1(phuang,y,pbis(i,j,k),'pchip');
            end
        end
    end
end

%% Compare spiciness 
figure
subplot(2,1,1)
pcolor(latitude,-zlev,squeeze(tau_remi(661,:,:))');shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(a) QL Pot. Spiciness - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
set(gca,'FontSize',10)
%pbaspect([2.4 1 1])
hold on
[C5,h5] = contour(latitude,-zlev,squeeze(tau_remi(661,:,:))', ...
    -1:0.1:1,'k', 'ShowText','on')
clabel(C5,h5,'FontSize',8) 
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
subplot(2,1,2)
pcolor(latitude,-zlev,squeeze(tau_remi_ref(661,:,:))');shading flat;
colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(b) QL Ref. Pot. Spiciness - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
set(gca,'FontSize',10)
%pbaspect([2.4 1 1])
hold on 
[C4,h4] = contour(latitude,-zlev,squeeze(tau_remi_ref(661,:,:))',  ...
    -1:0.1:1,'k', 'ShowText','on')
clabel(C4,h4,'FontSize',8) 
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
print('fig08_10pt.eps','-depsc2')

%% Compare patched potential spiciness versus refernece pot spicinses
figure
subplot(2,1,1)
pcolor(latitude,-zlev,squeeze(tau_ppd(661,:,:))');shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(a) Patched Pot. Spiciness - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
set(gca,'FontSize',10)
%pbaspect([2.4 1 1])
hold on
[C5,h5] = contour(latitude,-zlev,squeeze(tau_ppd(661,:,:))', ...
    -1:0.1:1,'k', 'ShowText','on')
clabel(C5,h5,'FontSize',8) 
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
subplot(2,1,2)
pcolor(latitude,-zlev,squeeze(tau_remi_ref(661,:,:))');shading flat;
colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(b) QL Ref. Pot. Spiciness - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
set(gca,'FontSize',10)
%pbaspect([2.4 1 1])
hold on 
[C4,h4] = contour(latitude,-zlev,squeeze(tau_remi_ref(661,:,:))',  ...
    -1:0.1:1,'k', 'ShowText','on')
clabel(C4,h4,'FontSize',8) 
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
print('fig08_10pt.eps','-depsc2')

%% Compare spicinesses
figure
subplot(3,1,1)
pcolor(latitude,-zlev,squeeze(tau_jmd(661,:,:))');shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(a) JMD Pot. Spiciness - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
hold on
contour(latitude,-zlev,squeeze(tau_jmd(661,:,:))',-1:0.1:1,'k', ...
    'ShowText','on')
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
subplot(3,1,2) 
pcolor(latitude,-zlev,squeeze(tau_remi(661,:,:))');shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(b) QL Pot. Spiciness - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
hold on
contour(latitude,-zlev,squeeze(tau_remi(661,:,:))',-1:0.1:1,'k', ...
    'ShowText','on')
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
subplot(3,1,3)
pcolor(latitude,-zlev,squeeze(tau_remi_ref(661,:,:))');shading flat;
colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(c) QL Ref. Pot. Spiciness - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
hold on 
contour(latitude,-zlev,squeeze(tau_remi_ref(661,:,:))',-1:0.1:1,'k', ...
    'ShowText','on')
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
subplot(3,1,3) 
print('fig08_bis.eps','-depsc2')

%% Compare potential spicities
figure
subplot(2,1,1)
pcolor(latitude,-zlev,squeeze(tau_huang(661,:,:))');shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(a) Pot. Spicity (p_r = 1000 dbar) - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
hold on
[C6,h6] = contour(latitude,-zlev,squeeze(tau_huang(661,:,:))', ...
    -1:0.1:1,'k','ShowText','on')
clabel(C6,h6,'FontSize',8)
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
subplot(2,1,2)
pcolor(latitude,-zlev,squeeze(pi_ref(661,:,:))');shading flat;
colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(b) Ref. Pot. Spicity - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
hold on 
[C7,h7]= contour(latitude,-zlev,squeeze(pi_ref(661,:,:))',-1:0.1:1,'k', ...
    'ShowText','on')
clabel(C7,h7,'FontSize',8) 
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
print('fig05_10pt.eps','-depsc2')

%% Compare potential ref spicity with patched potential spicity 
figure
subplot(2,1,1)
pcolor(latitude,-zlev,squeeze(pi_ppd(661,:,:))');shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(a) Patched Pot. Spicity - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
hold on
[C6,h6] = contour(latitude,-zlev,squeeze(pi_ppd(661,:,:))', ...
    -1:0.1:1,'k','ShowText','on')
clabel(C6,h6,'FontSize',8)
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
subplot(2,1,2)
pcolor(latitude,-zlev,squeeze(pi_ref(661,:,:))');shading flat;
colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['(b) Ref. Pot. Spicity - 30' char(176) ' W']) 
cmocean('balance')
caxis([-3 3])
hold on 
[C7,h7]= contour(latitude,-zlev,squeeze(pi_ref(661,:,:))',-1:0.1:1,'k', ...
    'ShowText','on')
clabel(C7,h7,'FontSize',8) 
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
axis([-80 70 -6000 0])
hold off
print('fig05_10pt.eps','-depsc2')

%% Potential spicity referenced 
figure
pcolor(latitude,-zlev,squeeze(tau_huang(661,:,:))');shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)')
title(['Potential Spicity (p_r = 1000 dbar) - 30' char(176) ' W'])
cmocean('balance')
caxis([-3 3])
hold on
contour(latitude,-zlev,squeeze(tau_huang(661,:,:))',-1:0.1:1,'k', ...
    'ShowText','on')
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

%% Make plot of potential spiciness along with sigma_1 contours
%plot_figure_introduction
plot_figure_introduction_revised

%% Make plots of the various spiciness/density correlations
plot_figure_correlations

%% valid points along 330W section
valid661 = (lon3d == 330.) & valid; 

%% Construction of potential spiciness anomaly
tau_anomaly = tau_remi_ref;
tau_anomaly(valid) = tau_remi_ref(valid) - ftau(gammat(valid)); 
rms_tau = rms(tau_anomaly(valid661));
tau_anomaly(valid) = tau_anomaly(valid)/rms_tau; 
figure
pcolor(latitude,-zlev,squeeze(tau_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(a) \tau_R anomaly - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2. 2.])
hold on
contour(latitude,-zlev,squeeze(tau_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on') 
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('tau_anomaly_gammat.eps','-depsc2')

%% Construction of spicity anomaly
pi_anomaly = tau_huang;
pi_anomaly(valid) = tau_huang(valid) - fpi(gammat(valid));
rms_pi = rms(pi_anomaly(valid661)); 
pi_anomaly(valid) = pi_anomaly(valid)/rms_pi;
figure
pcolor(latitude,-zlev,squeeze(pi_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel('Latitude (depth)');ylabel('Depth (m)') 
title('(b) Potential spicity anomaly - 30W')
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2 2])
hold on
contour(latitude,-zlev,squeeze(pi_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on') 
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('pi_anomaly_gammat.eps','-depsc2')

%% New construction of spicity anomaly
pi_anomaly = pi_ref;
pi_anomaly(valid) = pi_ref(valid) - fpi(gammat(valid));
rms_pi = rms(pi_anomaly(valid661)); 
pi_anomaly(valid) = pi_anomaly(valid)/rms_pi;
figure
pcolor(latitude,-zlev,squeeze(pi_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(b) \pi_R anomaly - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2 2])
hold on
contour(latitude,-zlev,squeeze(pi_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on') 
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('pi_anomaly_gammat.eps','-depsc2')

%% Construction of CT anomaly
ct_anomaly = ct;
ct_anomaly(valid) = ct(valid) - fct(gammat(valid));
rms_ct = rms(ct_anomaly(valid661)); 
ct_anomaly(valid) = ct_anomaly(valid)/rms_ct;
figure
pcolor(latitude(1:280),-zlev,squeeze(ct_anomaly(661,1:280,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(d) \Theta anomaly - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2. 2.])
hold on
contour(latitude,-zlev,squeeze(ct_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on') 
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('ct_anomaly_gammat.eps','-depsc2')

%% COnstruction of SR anomaly
sr_anomaly = sr;
sr_anomaly(valid) = sr(valid) - fsr(gammat(valid));
rms_sr = rms(sr_anomaly(valid661));
sr_anomaly(valid) = sr_anomaly(valid)/rms_sr;
figure
pcolor(latitude,-zlev,squeeze(sr_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(c) S_A anomaly - 30' char(176) ' W'])
axis([-80 70 -6000 0]) 
cmocean('balance') 
caxis([-2. 2.])
hold on
contour(latitude,-zlev,squeeze(sr_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on') 
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('sr_anomaly_gammat.eps','-depsc2') 

%% Combine anomaly plots
figure
subplot(2,2,1)
pcolor(latitude,-zlev,squeeze(tau_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(a) \tau''_{ref}/\Delta \tau - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2. 2.])
hold on
[C5,h5]=contour(latitude,-zlev,squeeze(tau_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
clabel(C5,h5,'FontSize',8)
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

subplot(2,2,2)
pcolor(latitude,-zlev,squeeze(pi_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(b) \pi''_{ref}/\Delta \pi - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2 2])
hold on
[C5,h5]= contour(latitude,-zlev,squeeze(pi_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
clabel(C5,h5,'FontSize',8) 
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

subplot(2,2,3)
pcolor(latitude,-zlev,squeeze(sr_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(c) S''_A/\Delta S - 30' char(176) ' W'])
axis([-80 70 -6000 0]) 
cmocean('balance') 
caxis([-2. 2.])
hold on
[C5,h5] = contour(latitude,-zlev,squeeze(sr_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
clabel(C5,h5,'FontSize',8)
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
set(gca,'FontSize',10) 
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off

subplot(2,2,4)
pcolor(latitude(1:280),-zlev,squeeze(ct_anomaly(661,1:280,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(d) \Theta''/\Delta \Theta - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2. 2.])
hold on
[C5,h5]= contour(latitude,-zlev,squeeze(ct_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
clabel(C5,h5,'FontSize',8)
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('fig12_10pt','-depsc2') 

%% Individual panels of Figure 12
figure
pcolor(latitude,-zlev,squeeze(tau_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(a) \tau''_{ref}/\Delta \tau - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2. 2.])
hold on
[C5,h5]=contour(latitude,-zlev,squeeze(tau_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
%clabel(C5,h5,'FontSize',8)
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
%set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('fig12panelA','-depsc2')

%%
figure
pcolor(latitude,-zlev,squeeze(pi_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(b) \pi''_{ref}/\Delta \pi - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2 2])
hold on
[C5,h5]= contour(latitude,-zlev,squeeze(pi_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
%clabel(C5,h5,'FontSize',8) 
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
%set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('fig12panelB','-depsc2')

%%
figure
pcolor(latitude,-zlev,squeeze(sr_anomaly(661,:,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(c) S''_A/\Delta S - 30' char(176) ' W'])
axis([-80 70 -6000 0]) 
cmocean('balance') 
caxis([-2. 2.])
hold on
[C5,h5] = contour(latitude,-zlev,squeeze(sr_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
%clabel(C5,h5,'FontSize',8)
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
%set(gca,'FontSize',10) 
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('fig12panelC','-depsc2')

%%
figure
pcolor(latitude(1:280),-zlev,squeeze(ct_anomaly(661,1:280,:))'); ...
    shading flat;colorbar
xlabel(['Latitude (' char(176) ' N)']);ylabel('Depth (m)') 
title(['(d) \Theta''/\Delta \Theta - 30' char(176) ' W'])
axis([-80 70 -6000 0])
cmocean('balance') 
caxis([-2. 2.])
hold on
[C5,h5]= contour(latitude,-zlev,squeeze(ct_anomaly(661,:,:))', ...
    -4:0.1:4,'k','ShowText','on');
%clabel(C5,h5,'FontSize',8)
contour(latitude,-zlev,squeeze(gammat(661,:,:))', ...
    27:0.1:30,'Color',[0.6350 0.0780 0.1840])
%set(gca,'FontSize',10)
fill(latlat,ztopzlow,[0.3 0.3 0.3])
hold off
print('fig12panelD','-depsc2')



%% Compute correlations
xoo = -4:0.1:4;
figure
subplot(2,2,1)
plot(squeeze(sr_anomaly(661,1:280,:)),squeeze(ct_anomaly(661,1:280,:)),'b.')
hold on
plot(xoo,xoo,'r')
grid on
hold off
axis([-4 4 -4 4])
xlabel('S''_A/\Delta S');ylabel('\Theta''/\Delta \Theta')
title('(a)')
subplot(2,2,2)
plot(squeeze(sr_anomaly(661,1:280,:)),squeeze(pi_anomaly(661,1:280,:)),'b.')
hold on
plot(xoo,xoo,'r')
grid on
hold off
axis([-4 4 -4 4])
title('(b)')
xlabel('S''_A/\Delta S');ylabel('\pi''_{ref}/\Delta \pi')
subplot(2,2,3)
plot(squeeze(sr_anomaly(661,1:280,:)),squeeze(tau_anomaly(661,1:280,:)),'b.')
hold on
plot(xoo,xoo,'r')
grid on
hold off
axis([-4 4 -4 4])
title('(c)') 
xlabel('S''_A/\Delta S');ylabel('\tau''_{ref}/\Delta \tau')
subplot(2,2,4)
plot(squeeze(ct_anomaly(661,1:280,:)),squeeze(pi_anomaly(661,1:280,:)),'b.')
hold on
plot(xoo,xoo,'r')
grid on
hold off
axis([-4 4 -4 4])
title('(d)')
xlabel('\Theta''/\Delta \Theta');ylabel('\pi''_{ref}/\Delta \pi')
print('fig13','-depsc2') 
print('fig13','-dpng')


%% Compute orthogonality 

grad_ct = zeros(npi,npj,npk,3);
[grad_ct(:,:,:,1),grad_ct(:,:,:,2),grad_ct(:,:,:,3)] = ...
    grad_density(ct,longitude,latitude,zlev,ra);
grad_ctan = zeros(npi,npj,npk,3);
[grad_ctan(:,:,:,1),grad_ctan(:,:,:,2),grad_ctan(:,:,:,3)] = ...
    grad_density(ct_anomaly,longitude,latitude,zlev,ra);
sin_ct = sin_angle(gradn_gt,grad_ct);
sin_ctan = sin_angle(gradn_gt,grad_ctan); 

grad_sr = zeros(npi,npj,npk,3);
[grad_sr(:,:,:,1),grad_sr(:,:,:,2),grad_sr(:,:,:,3)] = ...
    grad_density(sr,longitude,latitude,zlev,ra);
grad_sran = zeros(npi,npj,npk,3);
[grad_sran(:,:,:,1),grad_sran(:,:,:,2),grad_sran(:,:,:,3)] = ...
    grad_density(sr_anomaly,longitude,latitude,zlev,ra);
sin_sr = sin_angle(gradn_gt,grad_sr);
sin_sran = sin_angle(gradn_gt,grad_sran); 

grad_pi = zeros(npi,npj,npk,3);
[grad_pi(:,:,:,1),grad_pi(:,:,:,2),grad_pi(:,:,:,3)] = ...
    grad_density(pi_ref,longitude,latitude,zlev,ra);
grad_pian = zeros(npi,npj,npk,3);
[grad_pian(:,:,:,1),grad_pian(:,:,:,2),grad_pian(:,:,:,3)] = ...
    grad_density(pi_anomaly,longitude,latitude,zlev,ra);
sin_pi = sin_angle(gradn_gt,grad_pi);
sin_pian = sin_angle(gradn_gt,grad_pian); 


grad_tau = zeros(npi,npj,npk,3);
[grad_tau(:,:,:,1),grad_tau(:,:,:,2),grad_tau(:,:,:,3)] = ...
    grad_density(tau_remi_ref,longitude,latitude,zlev,ra);
grad_tauan = zeros(npi,npj,npk,3);
[grad_tauan(:,:,:,1),grad_tauan(:,:,:,2),grad_tauan(:,:,:,3)] = ...
    grad_density(tau_anomaly,longitude,latitude,zlev,ra);
sin_tau = sin_angle(gradn_gt,grad_tau);
sin_tauan = sin_angle(gradn_gt,grad_tauan); 

%% Mean properties and bar plots
mct = median(asin(sin_ct(:)),'omitnan')*180/pi;
mctan = median(asin(sin_ctan(:)),'omitnan')*180/pi;
msr = median(asin(sin_sr(:)),'omitnan')*180/pi;
msran = median(asin(sin_sran(:)),'omitnan')*180/pi;
mpi = median(asin(sin_pi(:)),'omitnan')*180/pi;
mpian = median(asin(sin_pian(:)),'omitnan')*180/pi;
mtau = median(asin(sin_tau(:)),'omitnan')*180/pi;
mtauan = median(asin(sin_tauan(:)),'omitnan')*180/pi;

mct661 = median(asin(sin_ct(661,:,:)),[2,3],'omitnan')*180/pi;
mctan661 = median(asin(sin_ctan(661,:,:)),[2,3],'omitnan')*180/pi;
msr661 = median(asin(sin_sr(661,:,:)),[2,3],'omitnan')*180/pi;
msran661 = median(asin(sin_sran(661,:,:)),[2,3],'omitnan')*180/pi;
mpi661 = median(asin(sin_pi(661,:,:)),[2 3],'omitnan')*180/pi;
mpian661 = median(asin(sin_pian(661,:,:)),[2 3],'omitnan')*180/pi;
mtau661 = median(asin(sin_tau(661,:,:)),[2 3],'omitnan')*180/pi;
mtauan661 = median(asin(sin_tauan(661,:,:)),[2 3],'omitnan')*180/pi;

subplot(2,1,1)
x = categorical({'S_A','\tau_{ref}','\pi_{ref}','\Theta'});
x = reordercats(x,{'S_A','\tau_{ref}','\pi_{ref}','\Theta'}); 
y = [msr msran; mtau mtauan; mpi mpian; mct mctan]; 
bar(x,y)
ylabel('Median angle (degrees)')
title('(a) Global Ocean')

subplot(2,1,2)
x = categorical({'S_A','\tau_{ref}','\pi_{ref}','\Theta'});
x = reordercats(x,{'S_A','\tau_{ref}','\pi_{ref}','\Theta'}); 
y = [msr661 msran661; mtau661 mtauan661; mpi661 mpian661; mct661 mctan661];
bar(x,y)
ylabel('Median angle (degrees)') 
title(['(b) Atlantic Ocean 30' char(176) ' W'])
print('fig02_bis','-depsc2') 


%% Redefine arrays 
[srg,ctg] = meshgrid(0:0.05:40,-2:0.05:40);
[zrefg,prefg,sigma_refg,gammat_refg]= gsw_gammat_analytic_os_2021(srg,ctg); 

phuang = [0.,500.,1000.,2000.,3000.,4000.,5000.];
y = phuang; 
validpig = (srg>10.&srg<40.) & (ctg>-2.&ctg<40.);
[nct,nsr] = size(srg); 
pig = zeros(nct,nsr); 

%%
disp('Compute Huang spicity') 
for j=1:nsr
    disp(j) 
    for i=1:nct
        if validpig(i,j)
            for l=1:7
                y(l) = gsw_pspi(srg(i,j),ctg(i,j),phuang(l));
            end
            pig(i,j)=interp1(phuang,y,prefg(i,j),'pchip');
        end
    end
end


%% Construct plots of anomalies in theta/S space
large_gt = gammat_refg>28.5;
small_gt = gammat_refg<23.; 
ctg_anomaly = ctg;
ctg_anomaly(:) = ctg(:) - fct(gammat_refg(:));
ctg_anomaly(large_gt) = NaN; ctg(large_gt) = NaN; 
gammat_refg(large_gt) = NaN; 
ctg_anomaly(small_gt) = NaN; ctg(small_gt) = NaN;
gammat_refg(small_gt) = NaN; 

srg_anomaly = srg;
srg_anomaly(:) = srg(:) - fsr(gammat_refg(:));
srg(large_gt) = NaN; srg(small_gt) = NaN; 
srg_anomaly(large_gt) = NaN; srg_anomaly(small_gt) = NaN;

pig_anomaly = pig;
pig_anomaly(:) = pig(:) - fpi(gammat_refg(:)); 
pig(large_gt) = NaN; pig(small_gt) = NaN;
pig_anomaly(large_gt) = NaN; pig_anomaly(small_gt) = NaN;

taug = spice_tailleux(srg,ctg,prefg,35.,20.);
taug_anomaly = taug;
taug_anomaly(:) = taug(:) - ftau(gammat_refg(:));
taug(large_gt) = NaN; taug(small_gt) = NaN;
taug_anomaly(large_gt) = NaN; taug_anomaly(small_gt) = NaN; 

%% Combined plots
figure
taug = spice_tailleux(srg,ctg,prefg,35.,20.);
taug_anomaly = taug;
taug_anomaly(:) = taug(:) - ftau(gammat_refg(:));
taug(large_gt) = NaN; taug(small_gt) = NaN;
taug_anomaly(large_gt) = NaN; taug_anomaly(small_gt) = NaN; 
subplot(2,2,1) 
[C10,h10]=contour(srg,ctg,taug,-10:0.5:40,'m--');
clabel(C10,h10,'Color','m','FontSize',8)
xlabel('S_A (g/kg)');ylabel(['\Theta (' char(176) ' C)'])
title('(a) \tau_{ref} versus \tau''_{ref}') 
hold on
[C11,h11]=contour(srg,ctg,gammat_refg,20:0.5:35,'Color', ...
    [0.6350 0.0780 0.1840]);
%clabel(C11,h11,'Color',[0.6350 0.0780 0.1840])
%clabel(C11,h11,'FontSize',8) 
[C12,h12]=contour(srg,ctg,taug_anomaly,-30:1:30,'k','ShowText','on');
clabel(C12,h12,'FontSize',8)
axis([33 40 -2 30])
set(gca,'FontSize',10)
hold off

pig_anomaly = pig;
pig_anomaly(:) = pig(:) - fpi(gammat_refg(:)); 
pig(large_gt) = NaN; pig(small_gt) = NaN;
pig_anomaly(large_gt) = NaN; pig_anomaly(small_gt) = NaN; 
subplot(2,2,2) 
[C9,h9]=contour(srg,ctg,pig,-10:0.5:40,'m--');
clabel(C9,h9,'Color','m','FontSize',8)
xlabel('S_A (g/kg)');ylabel(['\Theta (' char(176) ' C)'])
title('(b) \pi_{ref} versus \pi''_{ref}') 
hold on
[C8,h8]=contour(srg,ctg,gammat_refg,20:0.5:35,'Color',...
    [0.6350 0.0780 0.1840]);
%clabel(C8,h8,'Color',[0.6350 0.0780 0.1840])
%clabel(C8,h8,'FontSize',8)
[C10,h10]=contour(srg,ctg,pig_anomaly,-30:1:30,'k','ShowText','on')
clabel(C10,h10,'FontSize',8)
axis([33 40 -2 30])
set(gca,'FontSize',10)
hold off

srg_anomaly = srg;
srg_anomaly(:) = srg(:) - fsr(gammat_refg(:));
srg(large_gt) = NaN; srg(small_gt) = NaN; 
srg_anomaly(large_gt) = NaN; srg_anomaly(small_gt) = NaN;
subplot(2,2,3) 
contour(srg,ctg,srg,0:0.5:40,'m--')
xlabel('S_A (g/kg)');ylabel(['\Theta (' char(176) ' C)'])
title('(c) S_A versus S''_A') 
hold on
[C7,h7]=contour(srg,ctg,gammat_refg,20:0.5:35,'Color', ...
    [0.6350 0.0780 0.1840]);
clabel(C7,h7,'Color',[0.6350 0.0780 0.1840],'FontSize',8) 
[C8,h8]=contour(srg,ctg,srg_anomaly,-30:0.5:30,'k','ShowText','on');
clabel(C8,h8,'FontSize',8)
axis([33 40 -2 30])
set(gca,'FontSize',10)
hold off

subplot(2,2,4)
large_gt = gammat_refg>28.5;
small_gt = gammat_refg<23.; 
ctg_anomaly = ctg;
ctg_anomaly(:) = ctg(:) - fct(gammat_refg(:));
ctg_anomaly(large_gt) = NaN; ctg(large_gt) = NaN; 
gammat_refg(large_gt) = NaN; 
ctg_anomaly(small_gt) = NaN; ctg(small_gt) = NaN;
gammat_refg(small_gt) = NaN; 
%contour(srg,ctg,ctg,-2:2:30,'k--','ShowText','on')
[C6,h6]=contour(srg,ctg,ctg,-2:3:30,'--m');
%clabel(C6,h6,'Color','m')
%clabel(C6,h6,'FontSize',8) 
xlabel('S_A (g/kg)');ylabel(['\Theta (' char(176) ' C)'])
title('(d) \Theta versus \Theta''') 
hold on
[C5,h5]= contour(srg,ctg,gammat_refg,20:0.5:35,'Color', ...
    [0.6350 0.0780 0.1840]);
clabel(C5,h5,'Color',[0.6350 0.0780 0.1840],'FontSize',8) 
[C7,h7]=contour(srg,ctg,ctg_anomaly,-30:1:30,'k','ShowText','on')
clabel(C7,h7,'FontSize',8)
axis([33 40 -2 30])
set(gca,'FontSize',10)
hold off

print('fig11_10pt','-depsc2') 













































