function comp_surfstats(IBAS,CNAM,YEAR1,YEARN)

DIRA=[IBAS filesep CNAM filesep 'atm/hist'];
DIRO=[IBAS filesep CNAM filesep 'ocn/hist'];
ODIR='Data';
OFIL=[ODIR '/mht_' CNAM sprintf('_%4.4i-%4.4i.mat',YEAR1,YEARN)];
ANNUALOCN=false;

% earth radius
R=6371e3;

mkdir(ODIR);

areaatm=ncread('areaatm.nc','areacella');
areaocn=ncread('../micom_hydrography/Data/grid.nc','parea');

w=[31 28 31 30 31 30 31 31 30 31 30 31];
w=ones(1,12);
w=w/sum(w);

w_yr=1/(YEARN-YEAR1+1);
%
mhflx_atm_2d_ann=zeros(144,96);
mhflx_ocn=zeros(166,3);

for year=YEAR1:YEARN
  yyyy=sprintf('%4.4i',year);
  
  for month=1:12
    mm=sprintf('%2.2i',month);


    % read atmospheric fields 
    FLNT=ncread(fname,'FLNT');
    FSNT=ncread(fname,'FSNT');

    FLNS=ncread(fname,'FLNS');
    FSNS=ncread(fname,'FSNS');
 
    LHFLX=ncread(fname,'LHFLX');
    SHFLX=ncread(fname,'SHFLX');

%    PREC=(ncread(fname,'PRECSC')+ncread(fname,'PRECSL'))*3.337e8;
    PREC=0;
  
    mhflx_atm_2d_ann=mhflx_atm_2d_ann+w_yr*w(month)*(FSNT-FLNT-FSNS+FLNS+LHFLX+SHFLX+PREC);
  end

  % read ocean fields
  mhflx_ocn=mhflx_ocn+w_yr*ncread([DIRO filesep CNAM '.micom.hy.' yyyy '.nc'],'mhflx');
end

save zz_test.mat mhflx_atm_2d_ann;

%
lat_atm=ncread([DIRA filesep CNAM '.cam2.h0.' yyyy '-' mm '.nc'],'lat');
gw=ncread([DIRA filesep CNAM '.cam2.h0.' yyyy '-' mm '.nc'],'gw');
lat_ocn=ncread([DIRO filesep CNAM '.micom.hy.' yyyy '.nc'],'lat');
 

nlat=96;nlon=144;
pi=3.14159;

fluxx_ann=sum(mhflx_atm_2d_ann,1);

% integrate from the north
mhflx_atm_ann=zeros(length(lat_atm)+0,1);  % notice here: change with/without snow term

aaa=areaatm(1,:);

% for calculation without snow terms
if 1
for j=nlat:-1:1
 mhflx_atm_ann(j)=-sum(fluxx_ann(j:nlat).*aaa(j:nlat));
end
lat_atmm=lat_atm;
end

% for calculation with snow terms
if 0
lat_atmm(1)=-90;
for j=(1:nlat-1)
 lat_atmm(j+1)=(lat_atm(j+1)+lat_atm(j))/2;
 mhflx_atm_ann(j+1)=-sum(fluxx_ann(1:j).*aaa(1:j));
end
mhflx_atm_ann(nlat+1)==-sum(fluxx_ann(1:nlat).*aaa(1:nlat));
lat_atmm(end+1)=90;
end



save(OFIL,'mhflx_atm_ann','lat_atmm','mhflx_ocn','lat_ocn');
