clear all;
close all;
load grid;
%%
onem=9806;%pressure for one meter: 9806kg/m/s^2 (pasca)

FILE_FLX='~/output/H04S5CFC/cfcflx/H04S5CFC_flxacc_sum_1948-2007.nc';
FILE_CFC='/Users/heyc/Micom/ExpID/H04S5CFC/H04S5CFC_trc_2007.09.nc';
FILE_DZ='/Users/heyc/Micom/ExpID/H04S5CFC/H04S5CFC_3d_2007.09.nc';
flx1=nc_varget(FILE_FLX,'flx11_acc');
flx1(find(flx1==0))=NaN;
cfc11=nc_varget(FILE_CFC,'cfc11');
dz=nc_varget(FILE_DZ,'dp');

dz=dz/onem;
% flx1(isnan(flx1))=0;
% cfc11(find(isnan(cfc11)))=0;
% dz(find(isnan(dz)))=0;
%%

unix('ls ~/ExpID/H04S5CFC/H04S5CFC_flx_????.??.nc >filesflx');
filesflx=importdata('filesflx');
len=length(filesflx);
flx2=nc_varget(filesflx{1},'flx11_acc');
flx2=zeros(size(flx2));
for i=1:len
       flx2=flx2+nc_varget(filesflx{i},'flx11_acc');
end

%%

% int=0;
% for k=1:35
%     int=int+sum(sum(squeeze(cfc11(k,:,:)).*parea.*squeeze(dz(k,:,:))));
% end
int=cfc11.*dz;
int=sum(int,1);
int=squeeze(int);

diff1=flx1-int;
diff2=flx2-int;

%%
load grads;

figure;
Hm_axes=axesm('MapProjection','eqdcylin');
pcolorm(plat,plon,flx1);shading flat;
setm(gca,'Frame','on','MeridianLabel','on',...
	'MLabelParallel','south','MLabelLocation',60,...
	'ParallelLabel','on','PLabelLocation',30,...
	'MapLatLimit',[-90 90],'MapLonLimit',[0 360],'Origin',[0 180 0]);
axis tight;
title('flx1');
caxis([-0.5 3]*1e-5);
colormap(cmap);
cbfit(2);

figure;
Hm_axes=axesm('MapProjection','eqdcylin');
pcolorm(plat,plon,flx2);shading flat;
setm(gca,'Frame','on','MeridianLabel','on',...
	'MLabelParallel','south','MLabelLocation',60,...
	'ParallelLabel','on','PLabelLocation',30,...
	'MapLatLimit',[-90 90],'MapLonLimit',[0 360],'Origin',[0 180 0]);
axis tight;
title('flx2');
caxis([-0.5 3]*1e-5);
colormap(cmap);
cbfit(2);

figure;
Hm_axes=axesm('MapProjection','eqdcylin');
pcolorm(plat,plon,int);shading flat;axis tight;
setm(gca,'Frame','on','MeridianLabel','on',...
	'MLabelParallel','south','MLabelLocation',60,...
	'ParallelLabel','on','PLabelLocation',30,...
	'MapLatLimit',[-90 90],'MapLonLimit',[0 360],'Origin',[0 180 0]);
axis tight;
title('int');
caxis([-0.5 3]*1e-5);
colormap(cmap);
cbfit(2);

figure
Hm_axes=axesm('MapProjection','eqdcylin');
pcolorm(plat,plon,diff1);shading flat;axis tight;
setm(gca,'Frame','on','MeridianLabel','on',...
	'MLabelParallel','south','MLabelLocation',60,...
	'ParallelLabel','on','PLabelLocation',30,...
	'MapLatLimit',[-90 90],'MapLonLimit',[0 360],'Origin',[0 180 0]);
axis tight;
title('diff1:flx-int');
caxis([-0.5 3]*1e-5);
colormap(cmap);
cbfit(2);

figure
Hm_axes=axesm('MapProjection','eqdcylin');
pcolorm(plat,plon,diff2);shading flat;axis tight;
setm(gca,'Frame','on','MeridianLabel','on',...
	'MLabelParallel','south','MLabelLocation',60,...
	'ParallelLabel','on','PLabelLocation',30,...
	'MapLatLimit',[-90 90],'MapLonLimit',[0 360],'Origin',[0 180 0]);
axis tight;
title('diff2:flx-int');
caxis([-0.5 3]*1e-5);
colormap(cmap);
cbfit(2);

int0=nansum(nansum(int.*parea));
flx10=nansum(nansum(flx1.*parea));
flx20=nansum(nansum(flx2.*parea));

d1=int0-flx10;
d2=int0-flx20;