function comp_atmstats(IBAS,CNAM,varargin)

% set path names, etc
IDIR=[IBAS filesep CNAM filesep 'atm/hist']; 
VARS={'TREFHT','FSNT','FSNTC','FLNT','FLNTC','FSNS','FSNSC','FLNS','FLNSC'};
ODIR='Data';
YEAR1OLD=varargin{1};
YEAR1=varargin{1+length(varargin)-2};
YEARN=varargin{2+length(varargin)-2};

% create output folder
mkdir(ODIR);

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

% compute hemispheric ice extent and volume/mean thickness 
mskflg=0;
for year=YEAR1:YEARN 
  yyyy=sprintf('%4.4i',year);
  disp([yyyy])
  for month=1:12
    mm=sprintf('%2.2i',month);
    for VNAME=VARS
      VNAME1=char(VNAME);
      disp([IDIR filesep CNAM '.cam.h0.' yyyy '-' mm '.nc'])
      % check if file exists; if not, set output to nan
      FNAME=[IDIR filesep CNAM '.cam.h0.' yyyy '-' mm '.nc'];
      fid=fopen(FNAME,'r'); 
      if fid < 0 
        atmstats.(VNAME1)(month,year-YEAR1OLD+1)=nan;  
        disp(['Cannot read file ' FNAME])
        exit
      else 
        fclose(fid);
        fld=ncread([IDIR filesep CNAM '.cam.h0.' yyyy '-' mm '.nc'],VNAME1);     
        if mskflg == 0
          GFIL='remap_grid_fv1.9x2.5_20111104.nc';
          area=ncread(GFIL,'grid_area');
          areatot=sum(area(:));
          mskflg=1;
        end 
        atmstats.(VNAME1)(month,year-YEAR1OLD+1)=sum(fld(:).*area(:))/areatot;
      end
    end   
  end
end

% save statistics 
save([ODIR filesep 'atmstats_' CNAM '_' sprintf('%4.4i',YEAR1OLD) '-' sprintf('%4.4i',YEARN) '.mat'],...
  'atmstats','areatot'); 


