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];
NMEMTEST=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;



% read sensitivity hindcasts 
YEARSTEST=1985:1:1985;
boxst=zeros(10,length(YEARSTEST));
boxtt=zeros(10,length(YEARSTEST));
bmntt=zeros(10,length(YEARSTEST))+1e20;
bmxtt=zeros(10,length(YEARSTEST));
for s=1:length(YEARSTEST)
  s
  for m=1:NMEMTEST
    load(sprintf('Data/TS2D_noresm1-cmip6_icetest1_%4.4i1015_mem%2.2i_%4.4i-%4.4i_%4.4i-%4.4i.mat',YEARSTEST(s),m,ZTOP,ZBOT,YEARSTEST(s)+1,YEARSTEST(s)+10));
    S=reshape(S,[],size(S,3));
    T=reshape(T,[],size(T,3));
    boxst(:,s)=boxst(:,s)+S(ind,:)'*weights/NMEMTEST;
    boxtt(:,s)=boxtt(:,s)+T(ind,:)'*weights/NMEMTEST;
    bmntt(:,s)=min(bmntt(:,s),T(ind,:)'*weights);
    bmxtt(:,s)=max(bmxtt(:,s),T(ind,:)'*weights);
  end
end


save(sprintf('Data/TSBOX_cmip6_%s_%i-%i_test.mat',TAG,ZTOP,ZBOT),'boxst','boxtt','bmntt','bmxtt','NMEMTEST','YEARSTEST','volumetot','areatot')

end
end

