function [flag]=fun_avergedata(fn_in,fn_out,varname,yr_st,yr_end,scale_factor,units)
% clc;clear;close all;fclose all;
% fn_in='BB4CMIP-BC-biomass_1850-2015.nc';
% fn_out='test.nc';
% varname='BC';
% yr_st=1850;
% yr_end=1851;

%%
flag=0;

if isfile(fn_out)
    delete(fn_out);
end
% info=ncinfo(fn_in);

time=ncread(fn_in,'time');
time_bnds=ncread(fn_in,'time_bnds');
nt=numel(time);
time_unit=ncreadatt(fn_in,'time','units');
time_st_string=time_unit(11:end);
time_st=datenum(time_st_string);
% time_st_str=datestr(time_st,'yyyy-mm-dd HH:MM:SS');
DateVector = datevec(time_st);
year1=DateVector(1)+time_bnds(1,1)/365;
year2=year1+nt/12;

if yr_st<year1
    error('yr_st out of range!');
end

if yr_end>year2
    error('yr_end out of range!');
end

ind_st=(yr_st-year1)*12+1;
ind_end=(yr_end-year1)*12+12;
ind_count=ind_end-ind_st+1;

nyr=yr_end-yr_st+1;

start=[1 1 ind_st];
count=[inf inf ind_count];
var=ncread(fn_in,varname,start,count);
var(var>5.0000e+19)=0;
[n1,n2,n3]=size(var);
var_reshape=reshape(var,n1,n2,12,nyr);
var_avg=nanmean(var_reshape,4);
FillValue=1.0e20;
var_final=var_avg.*scale_factor;
var_final(isnan(var_final)==1)=FillValue;

lon=ncread(fn_in,'longitude');
if lon(1)>lon(end)
    var_final=flipud(var_final);
    lon=flipud(lon);
    disp('flipped lon!');
end
lat=ncread(fn_in,'latitude');
if lat(1)>lat(end)
    var_final=fliplr(var_final);
    lat=flipud(lat);
    disp('flipped lat!');
end

%%
time_new=time_bnds(1,1:12)'.*24;
date_st_str=[sprintf('%d',yr_st),'-01','-01',' 00:00:00'];
date_st_num=datenum(date_st_str);
for ii=1:12
    date_str=[sprintf('%d',yr_st),'-',sprintf('%02d',ii),'-01',' 00:00:00'];
    time_new(ii)=(datenum(date_str)-date_st_num)*24;
end

nccreate(fn_out,varname,'Dimensions',{'lon',n1,'lat',n2,'time',inf},...
    'Datatype','single','FillValue',FillValue,'ChunkSize',[n1 n2 12],'DeflateLevel',4);
% ncwrite(fn_out,varname,single(var_avg));
ncwrite(fn_out,varname,single(var_final));
var_long_name=ncreadatt(fn_in,varname,'long_name');
var_std_name=ncreadatt(fn_in,varname,'standard_name');
% var_units=ncreadatt(fn_in,varname,'units');
ncwriteatt(fn_out,varname,'long_name',var_long_name);
ncwriteatt(fn_out,varname,'standard_name',var_std_name);
ncwriteatt(fn_out,varname,'units',units);
ncwriteatt(fn_out,varname,'add_offset',0.0);
ncwriteatt(fn_out,varname,'scale_factor',1.0);
ncwriteatt(fn_out,varname,'missing_value',FillValue);
%%
nccreate(fn_out,'time','Dimensions',{'time',inf},...
    'Datatype','single','ChunkSize',12,'DeflateLevel',4);
ncwrite(fn_out,'time',single(time_new));
ncwriteatt(fn_out,'time','long_name','Time');
ncwriteatt(fn_out,'time','standard_name','Time');
ncwriteatt(fn_out,'time','units',['hours since ',date_st_str]);
ncwriteatt(fn_out,'time','calendar','standard');
ncwriteatt(fn_out,'time','axis','T');
%%
nccreate(fn_out,'lon','Dimensions',{'lon',n1},...
    'Datatype','single','DeflateLevel',4);
ncwrite(fn_out,'lon',single(lon));
ncwriteatt(fn_out,'lon','long_name','Longitude');
ncwriteatt(fn_out,'lon','units','degrees_east');
ncwriteatt(fn_out,'lon','axis','X');

%%
nccreate(fn_out,'lat','Dimensions',{'lat',n2},...
    'Datatype','single','DeflateLevel',4);
ncwrite(fn_out,'lat',single(lat));
ncwriteatt(fn_out,'lat','long_name','Latitude');
ncwriteatt(fn_out,'lat','units','degrees_north');
ncwriteatt(fn_out,'lat','axis','Y');

%%
ncwriteatt(fn_out,'/','creation_date',datestr(now));
ncwriteatt(fn_out,'/','Title','Global fire emissions from BB4CMIP6, monthly average');
ncwriteatt(fn_out,'/','Start Year',sprintf('%d',yr_st));
ncwriteatt(fn_out,'/','End Year',sprintf('%d',yr_end));
ncwriteatt(fn_out,'/','References','Van Marle et al. (2017). Historic global biomass burning emissions for CMIP6 (BB4CMIP) based on merging satellite observations with proxies and fire models (1750-2015). doi: https://doi.org/10.5194/gmd-10-3329-2017');
ncwriteatt(fn_out,'/','Conventions','COARDS');
ncwriteatt(fn_out,'/','Contact','Pengfei Liu (pengfeiliu@fas.harvard.edu)');



% info_out=ncinfo(fn_out);
% image(var_avg(:,:,7)'*1e13);
flag=1;


