function comp_ocnstats_ann(IBAS,CNAM,CTAG,varargin)

% set path names, etc
if strcmp(IBAS,'/work/ibe062/norstore9133')
  IDIR=[IBAS filesep CNAM filesep 'ocn'];
elseif strcmp(IBAS,'/hexagon/work/fanf/NorCPM')
  IDIR=[IBAS];
else
  IDIR=[IBAS filesep CNAM filesep 'ocn/hist'];
%  IDIR=[IBAS filesep CNAM filesep 'run/'];
end

GFIL='Data/grid_mis3.nc';  % !!!! NOTICE here!!!!
%GFIL='Data/grid.nc';        % !!!! NOTICE here!!!!

ODIR='Data';
VARS={'templvl','salnlvl'};
YEAR1OLD=varargin{1};
YEAR1=varargin{1+length(varargin)-2};
YEARN=varargin{2+length(varargin)-2};

% create output folder
mkdir(ODIR);

% read grid information and compute weights
lon=ncread(GFIL,'plon');
lat=ncread(GFIL,'plat');
mask=ncread(GFIL,'pmask');
area=ncread(GFIL,'parea');area(isnan(area))=0;
depth=ncread(GFIL,'pdepth');depth(isnan(depth))=0;

% read old file if exists 
if YEAR1OLD<YEAR1 
 load([ODIR filesep 'ocnstats_' CNAM CTAG '_' sprintf('%4.4i',YEAR1OLD) '-' sprintf('%4.4i',YEAR1-1) '.mat'],'ocnstats'); 
end

% compute nino3 index
weightsflg=0;
for year=YEAR1:YEARN 
  yyyy=sprintf('%4.4i',year);
  disp([yyyy]) 
    for VNAME=VARS
      FNAME=[IDIR filesep CNAM '.micom.hy.' yyyy '.nc']
      fid=fopen(FNAME,'r');
      if fid < 0
        ocnstats.global.(char(VNAME))(year-YEAR1OLD+1)=nan;
        ocnstats.profile.(char(VNAME))(year-YEAR1OLD+1,:)=nan;
      else 
        fclose(fid);
        if weightsflg == 0 
            weightsflg=1;
            depthbnds=ncread(FNAME,'depth_bnds');
            kdm=size(depthbnds,2);
            for k=1:kdm
              cellvolume(:,:,k)=area.*(min(depthbnds(2,k),depth)-min(depthbnds(1,k),depth));
              weights2d(:,:,k)=cellvolume(:,:,k)/max(1e-20,sum(sum(cellvolume(:,:,k))));
            end
            weights3d=cellvolume/sum(cellvolume(:));
        end 
        fld=ncread(FNAME,char(VNAME));
        fld(isnan(fld))=0;
        ocnstats.global.(char(VNAME))(year-YEAR1OLD+1)=sum(fld(:).*weights3d(:));
        ocnstats.profile.(char(VNAME))(year-YEAR1OLD+1,:)=squeeze(sum(sum(fld.*weights2d)));    
      end 
    end
end

for VNAME=VARS
 clim=nanmean(ocnstats.profile.(char(VNAME)));
 for k=1:size(ocnstats.profile.(char(VNAME)),2)
    ocnstats.profile_anomaly.(char(VNAME))(1:size(ocnstats.profile.(char(VNAME)),1),k)=ocnstats.profile.(char(VNAME))(:,k)-clim(k);
 end    
end

% add cellvolume
ocnstats.cellvolume=cellvolume;

% save index
save([ODIR filesep 'ocnstats_' CNAM CTAG '_' sprintf('%4.4i',YEAR1OLD) '-' sprintf('%4.4i',YEARN) '.mat'],'ocnstats'); 

