
%
%addpath ~/matlab/tools/m_map1.4
addpath ~/matlab/mats
addpath ~/matlab/tools

%
% Set a structure with experiments' information:
expr(1)=struct('name','NBF1850_f19_tn11_test_mis3b', ...
               'name_disp','MIS3', ...
               'first_year',1601, ...
               'last_year',1700, ...
               'display_first_year',1, ...
               'filetype','hm');


%
fgrid='/work/shared/noresm/inputdata/ocn/micom/tnx1v1mis3/20160721/grid.nc';
pclon=ncread(fgrid,'pclon');pclon=permute(pclon,[3 1 2]);
pclat=ncread(fgrid,'pclat');pclat=permute(pclat,[3 1 2]);
pmask=ncread(fgrid,'pmask');
parea=ncread(fgrid,'parea');
pdepth=ncread(fgrid,'pdepth');%pdepth(pdepth==0)=nan;
[nx ny]=size(pdepth);
plon=ncread(fgrid,'plon');%plon(pdepth==0)=nan;
plat=ncread(fgrid,'plat');%plat(pdepth==0)=nan;
xx=plon;xx=xx(250:end,:); xx(xx<0)=xx(xx<0)+360;plon(250:end,:)=xx;

%
ff=figure;clf;
set(gcf,'renderer','zbuffer');
set(ff,'paperunits','inches','paperposition',[1,1,6,6],'position',[100 100 1600 1500]);
set(0, 'DefaultAxesFontName', 'AvantGarde');hold on

for nexp=1:length(expr)

  fyear=expr(nexp).first_year;
  lyear=expr(nexp).last_year;
  expid=expr(nexp).name;
  ftype=expr(nexp).filetype

  % read in MIS3 processed output
  iceconc=load(['../iceconc_' expid '_' num2str(fyear) '-' num2str(lyear) '.mat']);
  iceconc.hice(iceconc.hice>10)=10;
  hice_mar=squeeze(iceconc.hice(:,:,3));
  fice_mar=squeeze(iceconc.fice(:,:,3));
  hice_sep=squeeze(iceconc.hice(:,:,9));
  fice_sep=squeeze(iceconc.fice(:,:,9));
  
  cc=8;
  hice_mar(hice_mar>cc)=cc;

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

  % first to plot MIS3 March ice thickness during 1901-2000
  m_proj('Stereographic','lon',-40,'lat',57,'rad',35);
%  m_proj('miller Cylindrical','lat',[40 80],'lon',[-70 20]);  
  m_contourf(lon',lat',hice_mar','levellist',-cc:0.1:cc,'linestyle','none');

  caxis([0 cc]);
  cmap=cbrewer('seq', 'YlGnBu',cc*10);cmap(1,:)=[1 1 1];colormap(cmap);
  cc=colorbar('h');set(cc,'position',[0.3 0.08 0.4 0.02],'fontsize',8);

  % ------ overlap 15% sea ice concentration ------
  hold on;
  m_contour(lon',lat',fice_mar','levellist',15,'color',0.6*ones(1,3),'linewidth',2);
  hold on;
  m_contour(lon',lat',fice_sep','levellist',15,'color',0.6*ones(1,3),'linestyle','--','linewidth',1.5);
  hold on;
  m_contour(plon',plat',pdepth','levellist',1,'color',0*ones(1,3),'linewidth',1,'color','r');
  hold on;
  
  % ------ overlap the grid points where fw is released
  
  % 1) along the Norwegian coast
  X=[113 113 113 114 114 115 116 117 118 118 118 118 118 118 119];
  Y=313:327;
  ind=(Y-1)*360+X;
  pclon_re=reshape(pclon,4,[]);pclat_re=reshape(pclat,4,[]);
  %hh=m_patch2(pclon_re(:,ind),pclat_re(:,ind),'k');set(hh,'edgecolor','k','linewidth',1,'facealpha',0.0);
  hold on;

  % 2) in the Nordic Seas
  fice_marr=fice_mar;fice_marr(:,end+1)=fice_marr(:,end);
  [X1 Y1]=find(plon>-20 & plon<20 & plat>64.2 & plat<76 & fice_marr<15);
  ind=[];ind=(Y1-1)*360+X1;
  pclon_re=reshape(pclon,4,[]);pclat_re=reshape(pclat,4,[]);
%  hh=m_patch2(pclon_re(:,ind),pclat_re(:,ind),'w');set(hh,'edgecolor','k','linewidth',1,'facealpha',0.0);

  % 2b) in the Nordic Seas
  iceconc2=load('../iceconc_NBF1850_f19_tn11_test_mis3b_1601-1700.mat');
  fice_sep2=squeeze(iceconc2.fice(:,:,9));
  fice_sepp=fice_sep2;fice_sepp(:,end+1)=fice_sepp(:,end);
  [X1b Y1b]=find(plon>-20 & plon<20 & plat>63.7 & plat<77 & fice_sepp<15);
  ind=[];ind=(Y1b-1)*360+X1b;
  pclon_re=reshape(pclon,4,[]);pclat_re=reshape(pclat,4,[]);
%  hh=m_patch2(pclon_re(:,ind),pclat_re(:,ind),'w');set(hh,'edgecolor','b','linewidth',1,'facealpha',0.0);

  % 2c) in the Nordic Seas; same with 2b but with extended grids in the south
  iceconc2=load('../iceconc_NBF1850_f19_tn11_test_mis3b_1601-1700.mat');
  fice_sep2=squeeze(iceconc2.fice(:,:,9));
  fice_sepp=fice_sep2;fice_sepp(:,end+1)=fice_sepp(:,end);
  [X1c Y1c]=find(plon>-20 & plon<20 & plat>63.8 & plat<77 & fice_sepp<15);
  XXX=[108:115 99:114 100:114 104:113 104:113 104:113 105:114];
  YYY=[318*ones(1,8) 317*ones(1,16) 316*ones(1,15) 315*ones(1,10) 314*ones(1,10) 313*ones(1,10) 312*ones(1,10)];
  X1c=[X1c;XXX'];Y1c=[Y1c;YYY'];
  ind=[];ind=(Y1c-1)*360+X1c;
  pclon_re=reshape(pclon,4,[]);pclat_re=reshape(pclat,4,[]);
  hh=m_patch2(pclon_re(:,ind),pclat_re(:,ind),'w');set(hh,'edgecolor','b','linewidth',0.5,'facealpha',0.0);


  % 3) in the western North Atlantic 
  fice_marr=fice_mar;fice_marr(:,end+1)=fice_marr(:,end);
  [X2 Y2]=find(plon>-66 & plon<-36 & plat>48 & plat<66 & fice_marr<15);
  ind=[];ind=(Y2-1)*360+X2;
  pclon_re=reshape(pclon,4,[]);pclat_re=reshape(pclat,4,[]);
  %hh=m_patch2(pclon_re(:,ind),pclat_re(:,ind),'w');set(hh,'edgecolor','k','linewidth',1,'facealpha',0.0);

  % 4) in the region between 50 and 70N
  fice_marr=fice_mar;fice_marr(:,end+1)=fice_marr(:,end);
  [X3 Y3]=find(plon>-70 & plon<20 & plat>50 & plat<70 & pdepth>0);
  ind=[];ind=(Y3-1)*360+X3;
  pclon_re=reshape(pclon,4,[]);pclat_re=reshape(pclat,4,[]);
%  hh=m_patch2(pclon_re(:,ind),pclat_re(:,ind),'w');set(hh,'edgecolor','c','linewidth',1,'facealpha',0.0);

  % 4) in the region between 40 and 55N
  fice_marr=fice_mar;fice_marr(:,end+1)=fice_marr(:,end);
  [X4 Y4]=find(plon>-80 & plon<0 & plat>40 & plat<55 & pdepth>0);
  ind=[];ind=(Y4-1)*360+X4;
  pclon_re=reshape(pclon,4,[]);pclat_re=reshape(pclat,4,[]);
%  hh=m_patch2(pclon_re(:,ind),pclat_re(:,ind),'w');set(hh,'edgecolor','c','linewidth',1,'facealpha',0.0);

lww=1.5;ccc='k';
% region 50-70 N
  lon1=-70:20;lat1=70*ones(1,length(lon1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);
  lon1=-70:20;lat1=50*ones(1,length(lon1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);
  lat1=50:70;lon1=-70*ones(1,length(lat1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);
  lat1=50:70;lon1=20*ones(1,length(lat1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);
% region 40-55 N
  ccc='b';
  lon1=-80:0;lat1=55*ones(1,length(lon1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);
  lon1=-80:0;lat1=40*ones(1,length(lon1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);
  lat1=40:55;lon1=-80*ones(1,length(lat1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);
  lat1=40:55;lon1=0*ones(1,length(lat1));
  m_line(lon1,lat1,'linewi',lww,'color',ccc);



  %
  m_coast('patch',[.7 .7 .7],'linewidth',0.2);
  %m_grid('fontsize',4,'xticklabel',[],'yticklabel',[]);
  m_grid('fontsize',8);
  title(expr(nexp).name_disp,'fontsize',4);

end

%return
print('-dpng','-r300','mask_fwf_tnx1v1.png');



%
% ------------------------------------------------------------------------------------
%
if 0

% to get the fwf_mask.nc file 
% u actually get a 'weight' for the fwf

% 1) get the mask along the norwegian coast
  fwf_mask=zeros(nx,ny);
  for nn=1:length(X)
    fwf_mask(X(nn),Y(nn))=1.0;
  end
  fwf_mask=fwf_mask/sum(parea(:).*fwf_mask(:));
%  fwf_mask=fwf_mask';
%
  fout='mask_fwf_tnx1v1.nc';
  delete(fout); 
  nccreate(fout,'fwf_mask','Dimensions', {'x', nx, 'y', ny});
  ncwrite(fout,'fwf_mask',fwf_mask,[1 1]);




% 2) get the mask for the Nordic Seas
fwf_mask=zeros(nx,ny);
for nn=1:length(X1)
  fwf_mask(X1(nn),Y1(nn))=1.0;
end
fwf_mask=fwf_mask/sum(parea(:).*fwf_mask(:));
%
fout='mask_fwf_nordic.nc';
delete(fout);
nccreate(fout,'fwf_mask','Dimensions', {'x', nx, 'y', ny});
ncwrite(fout,'fwf_mask',fwf_mask,[1 1]);

% 2b) get the mask for the Nordic Seas
fwf_mask=zeros(nx,ny);
for nn=1:length(X1b)
  fwf_mask(X1b(nn),Y1b(nn))=1.0;
end
fwf_mask=fwf_mask/sum(parea(:).*fwf_mask(:));
%
fout='mask_fwf_nordic_sep.nc';
delete(fout);
nccreate(fout,'fwf_mask','Dimensions', {'x', nx, 'y', ny});
ncwrite(fout,'fwf_mask',fwf_mask,[1 1]);


% 3) get the mask for the western North Atlantic
fwf_mask=zeros(nx,ny);
for nn=1:length(X2)
  fwf_mask(X2(nn),Y2(nn))=1.0;
end
fwf_mask=fwf_mask/sum(parea(:).*fwf_mask(:));
%
fout='mask_fwf_na.nc';
delete(fout);
nccreate(fout,'fwf_mask','Dimensions', {'x', nx, 'y', ny});
ncwrite(fout,'fwf_mask',fwf_mask,[1 1]);


% 4) get the mask for the region between 50 and 70N
fwf_mask=zeros(nx,ny);
for nn=1:length(X3)
  fwf_mask(X3(nn),Y3(nn))=1.0;
end
fwf_mask=fwf_mask/sum(parea(:).*fwf_mask(:));
%
fout='mask_fwf_50to70N.nc';
delete(fout);
nccreate(fout,'fwf_mask','Dimensions', {'x', nx, 'y', ny});
ncwrite(fout,'fwf_mask',fwf_mask,[1 1]);




% 5) get the mask for the region between 40 and 55N
fwf_mask=zeros(nx,ny);
for nn=1:length(X4)
  fwf_mask(X4(nn),Y4(nn))=1.0;
end
fwf_mask=fwf_mask/sum(parea(:).*fwf_mask(:));
%
fout='mask_fwf_40to55N.nc';
delete(fout);
nccreate(fout,'fwf_mask','Dimensions', {'x', nx, 'y', ny});
ncwrite(fout,'fwf_mask',fwf_mask,[1 1]);




% 2c) get the mask for the Nordic Seas
fwf_mask=zeros(nx,ny);
for nn=1:length(X1c)
  fwf_mask(X1c(nn),Y1c(nn))=1.0;
end
fwf_mask=fwf_mask/sum(parea(:).*fwf_mask(:));
%
fout='mask_fwf_nordic_sepL.nc';
delete(fout);
nccreate(fout,'fwf_mask','Dimensions', {'x', nx, 'y', ny});
ncwrite(fout,'fwf_mask',fwf_mask,[1 1]);


end
