% Plot the maps of upper ocean heat content inthe Western Pacific Ocean
% 
% Luc Rainville, December 2013.


% SSH (weekely)
aviso = load('./ssh_2001_2011.mat','ssh*');

% maps (monthly)
load('./maps_heat_content.mat')

% float data (every profile)
load('./KE_floats.mat','floats')
floats_all.time = floats.time;
floats_all.lon = floats.lon;
floats_all.lat = floats.lat;
floats_all.hc = floats.Q_500m;  % estimate of the heat content in the upper 500 m
clear floats

%% PLOT ONE EXAMPLE of a map of thickness


YEAR = 2005;
MONTH = 04;

XL1 = [125 185];
YL1 = [16 45];
   
lon = XL1(1):1/3:XL1(2);
lat =(YL1(1):1/3:YL1(2))';

N = 5*4;
dc = 0.2;
cmap = jet(N+2);
cmap(1,:)=[1 1 1]*1;
cmap(end,:)=[1 1 1]*0.5;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
m = find(TIME>datenum(YEAR,MONTH,15),1);

index_ssh = find(aviso.ssh_time>=TIME(m)-15 & aviso.ssh_time<TIME(m)+15);
ssh = interp2(aviso.ssh_lon,aviso.ssh_lat, nanmean(aviso.ssh_ssha(:,:,index_ssh),3)+aviso.ssh_mean2, lon, lat);

index = find(floats_all.time>=TIME(m)-15 & floats_all.time<TIME(m)+15 & ...
floats_all.lon>XL1(1) & floats_all.lon<XL1(2) & floats_all.lat>YL1(1) & floats_all.lat<YL1(2) & ...
isfinite(floats_all.hc) );

map0 = map_hc0(:,:,m)/1d10; % based on SSH only
map = map_hc(:,:,m)/1d10;   % with the floats
      
ii = find(map>=0);
i0 = find(map<0 | isnan(map));


% colorscale for scatter plot
c = floats_all.hc(index)/1d10;
clim = [0 N*dc]+54;
C = round((c-clim(1))/diff(clim)*(N-1))+1;
C(C<1)=1;
C(C>N)=N;
C=C+1; % I added a color at the beginning...

% for plotting
map(map-54<dc/2)=dc/2+54;
map(isnan(map))=-dc/2-54;
map(map-54> (N*dc-dc/2))=(N*dc-dc/2)+54;

figure(1)
clf

paper(10,5);
colormap(cmap);

pcolor(LON, LAT, map), shading flat;
hold on
contourf(bottom.lon, bottom.lat, bottom.H+2000,[1 1]*2000)
for n=1:length(index)
  plot(floats_all.lon(index(n)),floats_all.lat(index(n)),'sk','markersize',10,'linewidth',1,'markerfacecolor',cmap(C(n),:));
end
contour(lon, lat, ssh,-2:0.25:2,'color',[1 1 1]*0.5,'linewidth',1);
contour(lon, lat, ssh,[0.75 0.75],'color',[1 1 1]*0.5,'linewidth',2);
hold off
caxis([-dc N*dc+dc]+54)
colorbar;
xlim(XL1)
ylim(YL1)
set(gca,'xtick',120:10:180)
set(gca,'ytick',15:5:45)
set(gca,'fontsize',14);
title(['\bf',datestr(TIME(m),'mmm yyyy')])


%%
figure(2)
clf

clear mask
mask = map_hc_E(:,:,m)>0.3;
b = interp2(bottom.lon, bottom.lat, double(bottom.H>0), LON, LAT);
mask(b>0)=0;
mask=double(mask);

paper(10,5);
colormap(cmap);

contourf(LON, LAT, mask,[1 1]*0.5,'k')
caxis([-1 1])
xlim(XL1)
ylim(YL1)
set(gca,'xtick',120:10:180)
set(gca,'ytick',15:5:45)
set(gca,'fontsize',14);
title(['\bf MASK - ',datestr(TIME(m),'mmm yyyy')])

