clear all 


BOXLATVEC={[50 60],[50 60],[10 50],[48 65],[48 65]};
BOXLONVEC={[-45 -10],[-65 10],[-80 -5],[-60 -15],[-65 10]};
%ZBOTVEC=[1000 2000 10000];
ZBOTVEC=[2000];
NMEM=10;

for r=1:numel(BOXLATVEC)
BOXLAT=BOXLATVEC{r};
BOXLON=BOXLONVEC{r};
for k=1:numel(ZBOTVEC)
ZBOT=ZBOTVEC(k);


% set parameters 
%BOXLON=[0 20];BOXLAT=[65 80]; % Norwegian Sea 
%BOXLON=[1 3];BOXLAT=[65.5 66.5]; % Station M 
%BOXLON=[-45 -10];BOXLAT=[50 60]; % SPNA from Robson 2012 
%BOXLON=[-65 10];BOXLAT=[50 60]; % SPNA extended
%BOXLON=[-80 -5];BOXLAT=[10 50]; % STNA extended
ZTOP=0;
%ZBOT=1000; 
%ZBOT=10000; 

%BOXLON=[-60 -15];BOXLAT=[48 65];ZBOT=1000;
%BOXLON=[-60 -15];BOXLAT=[48 65];ZBOT=2000;
%BOXLON=[-60 -15];BOXLAT=[48 65];ZBOT=10000;
%BOXLON=[-65  10];BOXLAT=[48 65];ZBOT=1000;
%BOXLON=[-65  10];BOXLAT=[48 65];ZBOT=2000;
%BOXLON=[-65  10];BOXLAT=[48 65];ZBOT=10000;


% create box tag e.g. 45W10E50N60N 
TAG=[]
if BOXLAT(1)>=0 
   TAG=cat(2,TAG,int2str(BOXLAT(1)),'N');
else 
   TAG=cat(2,TAG,int2str(-BOXLAT(1)),'S');
end
if BOXLAT(2)>=0 
   TAG=cat(2,TAG,int2str(BOXLAT(2)),'N');
else 
   TAG=cat(2,TAG,int2str(-BOXLAT(2)),'S');
end
if BOXLON(1)>=0 
   TAG=cat(2,TAG,int2str(BOXLON(1)),'E');
else 
   TAG=cat(2,TAG,int2str(-BOXLON(1)),'W');
end
if BOXLON(2)>=0 
   TAG=cat(2,TAG,int2str(BOXLON(2)),'E');
else 
   TAG=cat(2,TAG,int2str(-BOXLON(2)),'W');
end

% read grid info
GRIDFILE='Data/grid_gx1.nc';
lon=ncread(GRIDFILE,'plon');
lat=ncread(GRIDFILE,'plat');
area=ncread(GRIDFILE,'parea');
mask=ncread(GRIDFILE,'pmask');
ind=find(lon>=BOXLON(1)&lon<=BOXLON(2)&lat>=BOXLAT(1)&lat<=BOXLAT(2)&mask>0.5);
load Data/maxdepth_model.mat
DZ=min(ZBOT,maxdepth(ind))-min(ZTOP,maxdepth(ind));
areatot=sum(area(ind));
volumetot=sum(DZ.*area(ind));
weights=DZ.*area(ind)./volumetot;


% process EN4 obs 
if true
YEAR1O=1950;
YEARNO=2015; 
%../../../Epocasa/SPNA-paper/Diagnose/HeatSalt/Data/TS2D_EN4_0000-2000_1950-2015.mat
load(sprintf('../../../Epocasa/SPNA-paper/Diagnose/HeatSalt/Data/TS2D_EN4_%4.4i-%4.4i_%4.4i-%4.4i.mat',ZTOP,ZBOT,YEAR1O,YEARNO));
[lato,lono]=meshgrid(-83:89,[1:359 0]); 
masko=ones(size(lono));
masko(find(S(:,:,1)<0.001))=0;
if BOXLON(2)<0 
indo=find(lono>=BOXLON(1)+360&lono<=BOXLON(2)+360&lato>=BOXLAT(1)&lato<=BOXLAT(2)&masko>0.5);
else
indo=find(lono>=BOXLON(1)&lono<=BOXLON(2)&lato>=BOXLAT(1)&lato<=BOXLAT(2)&masko>0.5);
end
weightso=cos(lato(indo)*pi/180)/sum(cos(lato(indo)*pi/180));
S=reshape(S,[],size(S,3));
T=reshape(T,[],size(T,3));
boxso=S(indo,:)'*weightso;
boxto=T(indo,:)'*weightso;
save('Data/comp_box_SPNA_EN4.mat','boxso','boxto','YEARNO')
end 

% read free historical 
YEAR1F=1950;
YEARNF=2010;
boxsf=zeros(YEARNF-YEAR1F+1,1);
boxtf=zeros(YEARNF-YEAR1F+1,1);
bmntf=zeros(YEARNF-YEAR1F+1,1)+1e20;
bmxtf=zeros(YEARNF-YEAR1F+1,1);
count=0;
for m=1:30
  m
  load(sprintf('Data/TS2D_noresm1-cmip6_historical_19500101_mem%2.2i_%4.4i-%4.4i_%4.4i-%4.4i.mat',m,ZTOP,ZBOT,YEAR1F,YEARNF));
  S=reshape(S,[],size(S,3));
  T=reshape(T,[],size(T,3));
  if sum(isnan(S(ind,:)+T(ind,:)))==0
    count=count+1
    boxsf=boxsf+S(ind,:)'*weights;
    boxtf=boxtf+T(ind,:)'*weights;
    bmntf=min(bmntf,T(ind,:)'*weights);
    bmxtf=max(bmxtf,T(ind,:)'*weights);
  end
end
boxsf=boxsf/count;
boxtf=boxtf/count;

% read analysis
if true
YEAR1A=1950;
YEARNA=2010; 
boxsa=zeros(YEARNA-YEAR1A+1,1);
boxta=zeros(YEARNA-YEAR1A+1,1);
bmnta=zeros(YEARNA-YEAR1A+1,1)+1e20;
bmxta=zeros(YEARNA-YEAR1A+1,1);
count=0; 
for m=1:30
  m
  load(sprintf('Data/TS2D_noresm1-cmip6_analysis_19500115_mem%2.2i_%4.4i-%4.4i_%4.4i-%4.4i.mat',m,ZTOP,ZBOT,YEAR1A,YEARNA));
  S=reshape(S,[],size(S,3));
  T=reshape(T,[],size(T,3));
  if sum(sum(isnan(S(ind,:)+T(ind,:))))==0 
    count=count+1
    boxsa=boxsa+S(ind,:)'*weights;
    boxta=boxta+T(ind,:)'*weights;
    bmnta=min(bmnta,T(ind,:)'*weights);
    bmxta=max(bmxta,T(ind,:)'*weights);
  end
end
boxsa=boxsa/count;
boxta=boxta/count;
end 

% read hindcasts 
YEARS=1960:1:2000;
boxsh=zeros(10,length(YEARS));
boxth=zeros(10,length(YEARS));
bmnth=zeros(10,length(YEARS))+1e20;
bmxth=zeros(10,length(YEARS));
for s=1:length(YEARS)
  s
  for m=1:NMEM 
    load(sprintf('Data/TS2D_noresm1-cmip6_hindcast_%4.4i1015_mem%2.2i_%4.4i-%4.4i_%4.4i-%4.4i.mat',YEARS(s),m,ZTOP,ZBOT,YEARS(s)+1,YEARS(s)+10));
    S=reshape(S,[],size(S,3));
    T=reshape(T,[],size(T,3));
    boxsh(:,s)=boxsh(:,s)+S(ind,:)'*weights/NMEM;
    boxth(:,s)=boxth(:,s)+T(ind,:)'*weights/NMEM;
    bmnth(:,s)=min(bmnth(:,s),T(ind,:)'*weights);
    bmxth(:,s)=max(bmxth(:,s),T(ind,:)'*weights);
  end   
end

save(sprintf('Data/TSBOX_cmip6_%s_%i-%i.mat',TAG,ZTOP,ZBOT),'boxsf','boxtf','bmntf','bmxtf','boxso','boxto','boxsa','boxta','bmnta','bmxta','boxsh','boxth','bmnth','bmxth','NMEM','YEAR1A','YEARNA','YEAR1F','YEARNF','YEAR1O','YEARNO','YEARS','volumetot','areatot')

% plot 
if false
figure(1)
clf 
subplot(2,1,1)
bar(YEAR1O:YEARNO,boxto-mean(boxto));
set(gca,'xlim',[1950 2020])
subplot(2,1,2)
bar(YEAR1A:YEARNA,boxta-nanmean(boxta));
set(gca,'xlim',[1950 2020])
OFFSET=nanmean(boxta);
hold on 
for s=1:length(YEARS)
  yvec=cat(1,boxta(YEARS(s)-YEAR1A),boxth(:,s))-OFFSET;
  xvec=YEARS(s)-1:YEARS(s)+9;
  plot(xvec,yvec,'r','linewidth',0.8);
end  
hold off

figure(2)
clf 
subplot(2,1,1)
bar(YEAR1O:YEARNO,boxso-mean(boxso));
set(gca,'xlim',[1950 2020])
subplot(2,1,2)
bar(YEAR1A:YEARNA,boxsa-nanmean(boxsa));
set(gca,'xlim',[1950 2020])
end

end 
end


