%
%addpath ~/matlab/m_map
addpath ~/matlab/mats
addpath ~/matlab/tools


% read in PI lon/lat ulon/ulat etc
%iceconcPI=load('../Data/iceconc_NBF1850OC_f19_tn11_02_1001-2000.mat');
iceconcPI=load('../Data/iceconc_NBF20TROC_f19_tn11_01_1976-2005.mat');

lon=double(iceconcPI.lon);lat=double(iceconcPI.lat);
x=lon;x=x(250:end,:); x(x<0)=x(x<0)+360;lon(250:end,:)=x;

% read in PI output
fice_marPI=squeeze(iceconcPI.fice(:,:,3));
fice_sepPI=squeeze(iceconcPI.fice(:,:,9));

hice_marPI=squeeze(iceconcPI.hice(:,:,3));
hice_sepPI=squeeze(iceconcPI.hice(:,:,9));

uvel_marPI=squeeze(iceconcPI.uvel(:,:,3));
vvel_marPI=squeeze(iceconcPI.vvel(:,:,3));
uvel_sepPI=squeeze(iceconcPI.uvel(:,:,9));
vvel_sepPI=squeeze(iceconcPI.vvel(:,:,9));


% read in HadISST data 
FNAME='../Data/HadISST_ice.nc';
lonn=double(ncread(FNAME,'lon')');
latt=double(ncread(FNAME,'lat')');
fice=permute(ncread(FNAME,'sea_ice_area_fraction'),[3 2 1]);
[lonn,latt]=meshgrid(lonn,latt);

% read in HadISST data 1970-1990
%marclim=squeeze(nanmean(fice((100*12+3):12:(120*12+3),:,:)));
%sepclim=squeeze(nanmean(fice((100*12+9):12:(120*12+9),:,:)));

% read in HadISST data 1976-2005
marclim=squeeze(nanmean(fice((105*12+3):12:(134*12+3),:,:)));
sepclim=squeeze(nanmean(fice((105*12+9):12:(134*12+9),:,:)));


%
%-----------------------------------------------------------------------------------------------------------
%
width=0.51;height=0.4;
pos_sp=[0.06 0.5 width height;
    0.52 0.5 width height;
    0.06 0.00 width height;
    0.52 0.00 width height];


ff=figure;
%set(ff,'paperunits','inches','paperposition',[1,1,6.6,1.9],'position',[100 100 1600 500]);
set(0, 'DefaultAxesFontName', 'AvantGarde');hold on

for nexp=1:4

  %
  % ------------------  SUBPLOT !! -----------------
  %

  subplot(2,2,nexp,'position',pos_sp(nexp,:))

  if nexp==1 
    m_proj('Stereographic','lon',-30,'lat',90,'rad',46);
    hice=hice_marPI;
    fice=fice_marPI;
    fice_obs=marclim;
    fout='a) NH Mar';
  elseif nexp==2
    m_proj('Stereographic','lon',-30,'lat',90,'rad',46);
    hice=hice_sepPI;
    fice=fice_sepPI;
    fice_obs=sepclim;
    fout='b) NH Sep';
  elseif nexp==3
    m_proj('Stereographic','lon',0,'lat',-90,'rad',40);
    hice=hice_marPI;
    fice=fice_marPI;
    fice_obs=marclim;
    fout='c) SH Mar';
  elseif nexp==4
    m_proj('Stereographic','lon',0,'lat',-90,'rad',40);
    hice=hice_sepPI;
    fice=fice_sepPI;
    fice_obs=sepclim;
    fout='d) SH Sep';
  end

  fice=fice/100;
  %

  m_contourf(lon',lat',hice','levellist',0:0.2:5,'linestyle','none');

  caxis([0 5]);
  cmap=vivid;colormap(cmap);
  cc=colorbar('v');set(cc,'position',[0.55 0.3 0.02 0.4],'fontsize',8,'ytick',0:5,'yAxisLocation','left');ylabel(cc,'(m)','fontsize',8);%title(cc,'(m)','fontsize',8)



  % ------ overlap observed 15% sea ice concentration ------
  hold on;
  m_contour(lon',lat',fice','levellist',0.15,'color','r','linestyle','-','linewidth',0.2);
  m_contour(lonn',latt',fice_obs','levellist',0.15,'color','k','linestyle','-','linewidth',0.7);
 
  %
  m_coast('patch',[.7 .7 .7],'linewidth',0.2);
  m_grid('fontsize',4,'ytick',[],'xtick',[],'xticklabel',[],'yticklabel',[]);
  title(fout,'fontsize',8);


end

%print(gcf,'-dpng','-r600',['ice_thickness_NBF1850OC_f19_tn11_02_1001-2000.png']);
%print(gcf,'-depsc','-r600',['ice_thickness_NBF1850OC_f19_tn11_02_1001-2000.eps']);

print(gcf,'-dpng','-r600',['ice_thickness_NBF20TROC_f19_tn11_01_1976-2005_2x2.png']);
print(gcf,'-depsc','-r600',['ice_thickness_NBF20TROC_f19_tn11_01_1976-2005_2x2.eps']);
