clear all 

% global parameters
FAC=100; % convert m to cm
LW=3;
FS=16;
AR=0.65;
YLABEL='1-2000m temperature anomaly (K)';
TAG='48N65N60W15W';
ZMAX=2000;
VERSION='cmip6';
STARTYEARS=[1986:1:1994];
%FIELD='T';CLIM=[-1 1]*0.5;
FIELD='S';CLIM=[-1 1]*0.1;

% read uninitialised-noAssim (TS2D_NorESM1-ME_historicalExt_noAssim_ensave_0000-2000_1950-2020.mat)
YEAR1F=1950;
YEARNF=2010; 
PREFIX='noresm1-cmip6_historical_19500101';
f=getfield(load(sprintf('Data/TS2D_%s_ensave_0000-%i_%i-%i.mat',PREFIX,ZMAX,YEAR1F,YEARNF),FIELD),FIELD);

% read observations 
[lato,lono]=meshgrid(-83:89,[1:359 0]);
YEAR1O=1950;
YEARNO=2015;
PREFIX='EN4';
o=getfield(load(sprintf('Data/TS2D_%s_0000-%i_%i-%i.mat',PREFIX,ZMAX,YEAR1O,YEARNO),FIELD),FIELD);
lono=cat(1,lono(end/2:end,:)-360,lono(1:end/2-1,:));
o=cat(1,o(end/2:end,:,:),o(1:end/2-1,:,:));
o=o(80:230,65:165,:);
o5010=mean(o(:,:,(1950-YEAR1O+1):(2010-YEAR1O)),3);

% read analysis 
  YEAR1A=1950;
  YEARNA=2010;
  PREFIX='noresm1-cmip6_analysis_19500115';
  a=getfield(load(sprintf('Data/TS2D_%s_ensave_0000-%i_%i-%i.mat',PREFIX,ZMAX,YEAR1A,YEARNA),FIELD),FIELD);

% read hindcasts and compute composite average
nstart=numel(STARTYEARS);
nyears=10; 
hc=zeros(size(a,1),size(a,2),nyears+1);
oc=zeros(size(o,1),size(o,2),nyears+1);
for n=1:nstart
  YEAR1H=STARTYEARS(n);
  YEARNH=YEAR1H+nyears-1;
  PREFIX='noresm1-cmip6_hindcast_';
  h=getfield(load(sprintf('Data/TS2D_%s%i_ensave_0000-%i_%i-%i.mat',PREFIX,YEAR1H-1,ZMAX,YEAR1H,YEARNH),FIELD),FIELD);
  h=cat(3,a(:,:,YEAR1H-YEAR1A),h);
  hc=hc+(h-f(:,:,(YEAR1H-YEAR1F):(YEARNH-YEAR1F+1)))/nstart;
  oc=oc+(o(:,:,(YEAR1H-YEAR1O):(YEARNH-YEAR1O+1))-o5010)/nstart;
end

% read grid and basin information 
GRIDDIR='/projects/NS2345K/noresm/inputdata/ocn/micom/gx1v6/20101119';
MERTRAOCEANS=[GRIDDIR filesep 'mertraoceans.dat'];
MERTRAINDEX=[GRIDDIR filesep 'mertraindex.dat'];
GRIDPATH=[GRIDDIR filesep 'grid.nc'];
mask=ncread(GRIDPATH,'pmask');
fid=fopen(MERTRAOCEANS);
dims=fscanf(fid,'%i',2);
fgetl(fid);
oceans=reshape(fscanf(fid,'%1i'),dims(2),dims(1))';
fclose(fid);

% mask land and oceans outside Atlantic Ocean 
for n=1:(nyears+1)
  h1=hc(:,:,n); 
  h1(find(abs(oceans-2)>0.5))=0;
  h1(find(mask<0.5))=nan;
  hc(:,:,n)=h1;
end

% plot 
figure(1)
clf
set(gcf,'inverthardcopy','off','color',[1 1 1]);
orient tall
addpath ~/matlab/cmap
load BlWhRe.mat
cmap(1,:)=[0.5 0.5 0.5];
%
subplot(4,2,1)
data=flipud(hc(:,:,1)');
imagesc(data(1:250,[270:end 1:60]));
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
ylabel('lead 0 yr','fontsize',14,'fontweight','bold');
%title('Hindcast - Free');
%
subplot(4,2,2)
data=flipud(oc(:,:,1)');
imagesc(data);
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
%ylabel('lead-yr 0');
%title('EN4 analysis');
%
subplot(4,2,3)
data=flipud(hc(:,:,2)');
imagesc(data(1:250,[270:end 1:60]));
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
ylabel('lead 1 yr','fontsize',14,'fontweight','bold');
%
subplot(4,2,4)
data=flipud(oc(:,:,2)');
imagesc(data);
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
%ylabel('lead-yr 1');
%
subplot(4,2,5)
data=flipud(mean(hc(:,:,3:6),3)');
imagesc(data(1:250,[270:end 1:60]));
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
ylabel('lead 2-5 yr','fontsize',14,'fontweight','bold');
%
subplot(4,2,6)
data=flipud(mean(oc(:,:,3:6),3)');
imagesc(data);
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
%ylabel('lead-yr 2-5');
%
subplot(4,2,7)
data=flipud(nanmean(hc(:,:,7:11),3)');
imagesc(data(1:250,[270:end 1:60]));
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
ylabel('lead 6-10 yr','fontsize',14,'fontweight','bold');
%
subplot(4,2,8)
data=flipud(mean(oc(:,:,7:11),3)');
imagesc(data);
caxis(CLIM)
colormap(cmap)
pos=get(gca,'position');
set(gca,'xtick',[],'ytick',[],'position',[pos(1) pos(2) pos(3)*1.2 pos(4)*1.2]);
%ylabel('lead-yr 6-10');
PREFIX=sprintf('Atlantic%s_0-%im_s%4.4i-%4.4i_%s',FIELD,ZMAX,STARTYEARS(1),STARTYEARS(end),VERSION);
print(gcf,'-depsc2',[PREFIX '.eps']);
print(gcf,'-dpng','-r150',[PREFIX '.png']);
