function comp_TS2D(IBAS,CNAM,ZTOP,ZBOT,YEAR1,YEARN) 

% set path names, etc
IDIR=[IBAS filesep CNAM filesep 'ocn/hist'];

GFIL='Data/grid_gx1.nc';
bath=ncread(GFIL,'pdepth');
[idm,jdm]=size(bath);
T=zeros(idm,jdm,YEARN-YEAR1+1,'single');
S=zeros(idm,jdm,YEARN-YEAR1+1,'single');
mweights=[31 28 31 30 31 30 31 31 30 31 30 31];mweights=mweights/sum(mweights);

% compute amo index
for year=YEAR1:YEARN 
  for month=1:12 
    FNAME=[IDIR filesep CNAM '.micom.hm.' sprintf('%4.4i-%2.2i',year,month) '.nc'];
    disp(FNAME)
    fid=fopen(FNAME,'r');
    if fid < 0
      T(:,:,year-YEAR1+1)=nan;
      S(:,:,year-YEAR1+1)=nan;
    else 
      if year==YEAR1 & month==1  
        depth=ncread(FNAME,'depth');
        depth_bnds=ncread(FNAME,'depth_bnds');
        data=ncread(FNAME,'templvl');  
        [idm,jdm,kdm]=size(data);
        weights3d=zeros(idm,jdm,kdm);
        ztop2d=ones(idm,jdm)*ZTOP;
        zbot2d=min(bath,ZBOT); 
        for k=1:kdm 
          ind=find(depth_bnds(1,k)<zbot2d);
          weights2d=zeros(idm,jdm);
          if ~isempty(ind)
            weights2d(ind)=(min(zbot2d(ind),depth_bnds(2,k))-depth_bnds(1,k))./(zbot2d(ind)-ztop2d(ind));     
          end
          weights3d(:,:,k)=weights2d;
        end         
        weights3dsuminv=max(1./nansum(weights3d,3),1e-20);
        for k=1:kdm
           weights3d(:,:,k)=weights3d(:,:,k).*weights3dsuminv;
        end
      end 
      fclose(fid);
      T(:,:,year-YEAR1+1)=T(:,:,year-YEAR1+1)+mweights(month)*...
        squeeze(nansum(ncread(FNAME,'templvl').*weights3d,3)); 
      S(:,:,year-YEAR1+1)=S(:,:,year-YEAR1+1)+mweights(month)*...
        squeeze(nansum(ncread(FNAME,'salnlvl').*weights3d,3)); 
    end
  end
end

% save index
save(['Data/TS2D_' CNAM '_' sprintf('%4.4i-%4.4i_%4.4i-%4.4i',ZTOP,ZBOT,YEAR1,YEARN) '.mat'],'T','S'); 

