% MeltPondGRLplotv3.m - script to plot figure 3 for Kim et al
% submission (2018) showing internal heating profile, equation of state for
% density as function of temperature, snapshots of T(x,z) in fresh and
% salty ponds, and width/time averaged T(z) for fresh and salty ponds

% code has only been tested on original author's system,
% and we accept no responsibility for adverse consequences if used on other systems.
% Use at your own risk.

% loads postprocessed data files from folder GRLplotdata/ - may need to
% unzip

% author: Andrew Wells, Oxford
% contact andrew.wells@physics.ox.ac.uk for further information


% set font sizes
axfont=16;
textfont=20;

% set default font sizes

set(0, 'DefaultTextFontName','Helvetica')
set(0, 'DefaultAxesFontName','Helvetica')
set(0, 'DefaultAxesFontSize',axfont)
set(0, 'DefaultTextFontSize',axfont)
set(0,'DefaultLineLineWidth',2);
set(0,'DefaultLineMarkerSize',8)

% generate figure
fig=figure(1);
scalfac=16.4/19;


set(gcf,'Units', 'centimeters','Position', [3 3 19*scalfac 14.25*scalfac],'PaperPosition', [0 0 19*scalfac 14.25*scalfac],'PaperSize',[19*scalfac 14.25*scalfac]);


hold on;

%%%%%%%%%%%%%% FIRST GENERATE AXES FOR CORE PLOTS %%%%
%
%


% set up panels 1, 2, 3 = a, b, c
% panels 4, 5, 6, 7 = d, e, f, g



% Q(z) internal heating plots for figure (a) 
load('GRLplotdata/Qz.mat')

figsub1=subplot('Position',[0.07 0.75 0.2 0.2]);
hold on
% plot heating profile, and label axes
plot(Gshape,Y,'-r')
axis([0 16 0 1])

box on
text(-5,1.15,'(a) Heating Profile','FontSize',textfont,'FontName','Helvetica','FontWeight','bold')
text(-4,0.5,'$\displaystyle{\frac{z}{H}}$','Rotation',0,'Interpreter','latex','FontSize',textfont)


set(gca,...
    'XTick',0:5:15,...
'YTick',0:0.5:1,...
'YTickLabel',{'0.0', '0.5', '1.0'},...
'FontUnits','points',...
'FontWeight','normal',...
'FontSize',axfont,...
'FontName','Helvetica')
xlabel({'$Q/Q_0$'},...
'Rotation',0,...
'FontUnits','points',...
'interpreter','latex',...
'FontSize',textfont,...
'FontName','Helvetica')



% load density as a function of temperature for fresh and salty ponds rho(T); 

load('GRLplotdata/rhoT.mat')

% plot figures (b) (c), labelling axes
figsub2=subplot('Position',[0.425 0.75 0.2 0.2]);
hold on
plot(T2,rho2,'-r');
axis([0 4 999.88 1000.04])

box on
text(-2.5,1000.065,'(b) Freshwater Density ','FontSize',textfont,'FontName','Helvetica','FontWeight','bold')
text(-1.75,999.975,'$\rho$','Rotation',0,'Interpreter','latex','FontSize',textfont)

text(-2.5,999.925,'$(\mathrm{kg\,m}^{-3})$','Interpreter','latex','FontSize',textfont)
set(gca,...
    'XTick',0:1:4,...
'YTick',999.9:0.05:1000.0,...
'YTickLabel',{'999.90', '999.95', '1000.00'},...
'FontUnits','points',...
'FontWeight','normal',...
'FontSize',axfont,...
'FontName','Helvetica')
xlabel({'$T, \quad(^{\circ}\mathrm{C})$'},...
'Rotation',0,...
'FontUnits','points',...
'interpreter','latex',...
'FontSize',textfont,...
'FontName','Helvetica')

% panel (c)
figsub3=subplot('Position',[0.77 0.75 0.2 0.2]);
hold on
plot(T3,rho3,'-r');
axis([-1.1 1.5 1016.03 1016.07])

box on
text(-2.5,1016.077,'(c) 20 g kg^{-1} Saltwater Density','FontSize',axfont,'FontSize',textfont,'FontName','Helvetica','FontWeight','bold')
text(-2.2,1016.057,'$\rho$','Rotation',0,'Interpreter','latex','FontSize',textfont)

text(-2.5,1016.042,'$(\mathrm{kg\,m}^{-3})$','Interpreter','latex','FontSize',textfont)
set(gca,...
    'XTick',-1:1:1,...
    'XTickLabel',{'-1', '0', '1'},...
'YTick',1016.03:0.02:1016.07,...
'FontUnits','points',...
'FontWeight','normal',...
'FontSize',axfont,...
'FontName','Helvetica')
xlabel({'$T, \quad(^{\circ}\mathrm{C})$'},...
'Rotation',0,...
'FontUnits','points',...
'interpreter','latex',...
'FontSize',textfont,...
'FontName','Helvetica')


% fresh pond, T(z) plot

% load temperature profile data for fresh pond 
load('GRLplotdata/Tzfresh.mat')
jmax=size(zstoreGRL,2);
col=coolwarmcmap;

% loop over all heating rates
for jj=1:1:jmax
    
   figsub5=subplot('Position',[0.77 0.425 0.16 0.2]);
%figsub5;
 hold on
 
 % plot T(z) for each heating rate (postprocessed from simulation output), with colour set by heating rate
%figure(2); 
plot(TdegCstoreGRL(:,jj),zstoreGRL(:,jj),'-','Color',col(colvec(jj),:))
hold on



end

% 
cbarhandle=colorbar;
colormap(coolwarmcmap)

% format axes, label colour bar

box on
axis([0 0.5 0 1])
text(-0.22,1.125,'(e)','FontSize',textfont,'FontName','Helvetica','FontWeight','bold')
text(-0.22,0.5,'$\displaystyle{\frac{z}{H}}$','Rotation',0,'Interpreter','latex','FontSize',textfont)
text(0.725,0.5,'$\overline{Q} H$','Interpreter','latex','FontSize',textfont)
text(0.695,0.3,'$(\mathrm{W\,m}^{-2})$','Interpreter','latex','FontSize',textfont)
set(gca,...
    'YTick',0:0.5:1,...
    'YTickLabel',{'0.0', '0.5', '1.0'},...
'XTick',0:0.2:0.4,...
'FontUnits','points',...
'FontWeight','normal',...
'FontSize',axfont,...
'FontName','Helvetica')
xlabel({'$T, \quad(^{\circ}\mathrm{C})$'},...
'Rotation',0,...
'FontUnits','points',...
'interpreter','latex',...
'FontSize',textfont,...
'FontName','Helvetica')


% set colour bar and ticks
set(cbarhandle,'Ytick',(colvecQHc-1)/256)
set(cbarhandle,'YtickLabel',QHcbarlab)
set(cbarhandle,'FontName','Helvetica','FontSize',axfont)


% panel f/g

% load temperature data for the 20psu case
load('GRLplotdata/TforS20.mat')



% panel g - T(z) for different S
figsub7=subplot('Position',[0.77 0.1 0.2 0.2]);

% % plot convective and conductive solutions for saline pond


plot(TconvS20F0p5,zconvS20F0p5,'-k')
hold on
% plot conductive solutions
plot(TcondS20F0p25,zcondS20F0p25,'--','Color',[0.5 0.5 0.5])
plot(TcondS20F0p25,zcondS20F0p25,'--r')
plot(TcondS20F0p5,zcondS20F0p5,'--m')
plot(TcondS20F0p75,zcondS20F0p75,'--b')


% generate legend
leg7=legend('Convective','Conductive');


 set(leg7,'Position',[0.915 0.12 0.035 0.025])
 set(leg7,'FontSize',18,'FontName','Helvetica')


% configure labels/annotations

box on
axis([-1.1 0.3 0 1])
text(-1.45,1.15,'(g)','FontSize',textfont,'FontName','Helvetica','FontWeight','bold')
text(-1.45,0.5,'$\displaystyle{\frac{z}{H}}$','Rotation',0,'Interpreter','latex','FontSize',textfont)

text(-0.75,0.875,'$f=0.75$','Interpreter','latex','FontSize',axfont,'Color','b')
text(-0.25,0.875,'$f=0.5$','Interpreter','latex','FontSize',axfont,'Color','m')
text(-0.15,0.4,'$f=0.25$','Interpreter','latex','FontSize',axfont,'Color','r')

set(gca,...
    'YTick',0:0.5:1,...
    'YTickLabel',{'0.0', '0.5', '1.0'},...
'XTick',-1:0.5:0.0,...
 'XTickLabel',{'-1.0', '-0.5', '0.0'},...
'FontUnits','points',...
'FontWeight','normal',...
'FontSize',axfont,...
'FontName','Helvetica')
xlabel({'$T, \quad(^{\circ}\mathrm{C})$'},...
'Rotation',0,...
'FontUnits','points',...
'interpreter','latex',...
'FontSize',textfont,...
'FontName','Helvetica')



% panel f - plot T(x,z)
figsub6=subplot('Position',[0.07 0.1 0.6/1.05 0.2/1.05]);


% colour filled plot of temperature T(x,z) snapshot
imagesc(xmapS,zmapS,TmapS)

% configure plot and colorbar

axis equal

    cbarhandleST=colorbar;
 
     caxis([-1.1 -0.8])
colormap(coolwarmcmap)


% set colour bar and ticks
set(cbarhandleST,'Ytick',[-1.1 -1.0 -0.9 -0.8],'FontSize',axfont)
set(cbarhandleST,'YtickLabel',{'-1.1','-1.0','-0.9','-0.8'})

% configure plot style and annotations
box on

text(-0.4,-0.25,'(f) 20 g kg^{-1} Saltwater Pond','FontSize',axfont,'FontSize',textfont,'FontName','Helvetica','FontWeight','bold')
text(-0.4,0.5,'$\displaystyle{\frac{z}{H}}$','Rotation',0,'Interpreter','latex','FontSize',textfont)

text(4.6,0.5,'$T$','Interpreter','latex','FontSize',textfont)
text(4.55,0.7,'$(^{\circ}\mathrm{C})$','Interpreter','latex','FontSize',textfont)
set(gca,...
    'YTick',0:0.5:1,...
    'YTickLabel',{'1.0', '0.5', '0.0'},...
'XTick',0:1:4,...
'FontUnits','points',...
'FontWeight','normal',...
'FontSize',axfont,...
'FontName','Helvetica')
xlabel({'$x/H$'},...
'Rotation',0,...
'FontUnits','points',...
'interpreter','latex',...
'FontSize',textfont,...
'FontName','Helvetica')


% panel d - plot T(x,z) for Fresh pond
load('GRLplotdata/TforFresh.mat')


figsub4=subplot('Position',[0.07 0.425 0.6/1.05 0.2/1.05]);

% color-filled plot of T(x,z) snapshot
imagesc(xmapF(1:512),zmapF,TmapF(:,1:512))

% configure plot style/annotations
axis equal

    cbarhandleST=colorbar;
   
     caxis([0.0 0.17])
colormap(coolwarmcmap)


% set colour bar and ticks
set(cbarhandleST,'Ytick',[0.0 0.05 0.1 0.15],'FontSize',axfont)
set(cbarhandleST,'YtickLabel',{'0.00','0.05','0.10','0.15'})

box on

text(-0.4,-0.25,'(d) Freshwater Pond','FontSize',axfont,'FontSize',textfont,'FontName','Helvetica','FontWeight','bold')
text(-0.4,0.5,'$\displaystyle{\frac{z}{H}}$','Rotation',0,'Interpreter','latex','FontSize',textfont)

text(4.6,0.5,'$T$','Interpreter','latex','FontSize',textfont)
text(4.55,0.7,'$(^{\circ}\mathrm{C})$','Interpreter','latex','FontSize',textfont)
set(gca,...
    'YTick',0:0.5:1,...
    'YTickLabel',{'1.0', '0.5', '0.0'},...
'XTick',0:1:4,...
'FontUnits','points',...
'FontWeight','normal',...
'FontSize',axfont,...
'FontName','Helvetica')
xlabel({'$x/H$'},...
'Rotation',0,...
'FontUnits','points',...
'interpreter','latex',...
'FontSize',textfont,...
'FontName','Helvetica')




% Finally print figures plots
print(fig,'PlottingDir/Figure3.eps','-depsc')
print(fig,'PlottingDir/Figure3.pdf','-dpdf')

savefig(fig,'PlottingDir/Figure3.fig')
