%atn 16.jan.2021
%script to mimic eccov4r3 heat budget closure, using ASTEr1 binary/netcdf outputs
%this code has been verified to work for test3d=1, not yet tested for depth-integrated.
clear all;
warning off

%some numbers relevant for ASTE domain
nx=270;nz=50;ncut1=450;ncut2=180;
ny=2*ncut1+nx+ncut2;nfx=[nx 0 nx ncut2 ncut1];nfy=[ncut1 0 nx nx nx];
fac=270/nx; 

RunStr='OFFICIAL_ASTE_R1';

fprintf('NEED TO set up dirRoot and subdir prior to proceeding\n');
dirRoot='';%['./' RunStr '/'];
if(isempty(dirRoot));error('edit script to set up dirRoot and subdir');end;

%set up directory, here called dirRoot, where everything is put in here including:
% dirRoot/example_budget_tracers
% dirRoot/diags/[STATE,TRSP,BUDG]/ if downloading raw binary files
% dirRoot/[nctiles_monthly,nctiles_monthly_integrated,nctiles_monthly_snapshots,nctiles_monthly_snapshots_integrated]
%
dirCur=[dirRoot '/example_budget_tracers/'];cd(dirCur);
addpath([dirRoot '/ASTE_mlibrary/']);			%<-- define your own location
addpath([dirRoot '/ASTE_mlibrary/gcmfaces_mod/']);
dirOut=dirCur;

%choose to read in raw binary or netcdf
fprintf('change read_netcdf=0/1 to switch between read netcdf and binary, default is 1\n');
read_netcdf=1;

if(read_netcdf);dirGrid=[dirRoot '/nctiles_grid/'];else;dirGrid=[dirRoot '/GRID_real8/'];end;

if(read_netcdf);dirIn=[dirRoot];else;dirIn=[dirRoot '/diags/BUDG/'];end;
dirMat=[dirRoot '/mat/diags_set_D/'];if(exist(dirMat)==0);mkdir(dirMat);end;

strbudg='Heat';
kBudget=1;
test3d=1;
save_budg_scalar=0;save_budg_lev=0; %<-- these were generally not used
save_budg_3d=0;save_budg_2d=0;	    %<-- typically save budg_2d.  budg_3d is VERY LARGE
plot_fig=1;
kz=[1 5;6 10;11 19;20 29;30 37;38 50];	%pick some range of depth level for plotting

% Create myparms structure
myparms=[];
myparms.yearFirst=2002;
myparms.yearLast=2017;
myparms.yearInAv=[2002 2017];
myparms.timeStep=600;%600 for NLFS+RFWF, 1200 for LFS+virtual salt
myparms.iceModel=1;
myparms.useRFWF=1;%0;
myparms.useNLFS=4;%0;
myparms.rhoconst=1029;
myparms.rcp=4109826;
myparms.rhoi=910;
myparms.rhosn=330;
myparms.flami=334000;
myparms.flamb=2500000;
myparms.SIsal0=4;
myparms.diagsNbRec=3;
myparms.diagsAreMonthly=1;
myparms.diagsAreAnnual=0;
myparms.recInAve=[1 192];
myparms.rStar=2;%0;
myparms.SPHF=0;

myparms.conserveTr=0;

deltaTime=myparms.timeStep;

%get time-steps:
if(read_netcdf);
  flist=load([dirIn 'time_list.txt']);	%<-- this file is in doc/ and also inside each nctiles_* dir
else;
  flist=dir([dirIn 'budg2d_snap_set2.*.data']);
  idot=find(flist(1).name=='.');idot=idot(1)+1:idot(2)-1;
end;

%set up mygrid structure in native binary format:
							%/*{{{*/ [grid]
global mygrid
mygrid=[];
mygrid.dirGrid=dirGrid;
mygrid.nFaces=5;
mygrid.fileFormat='compact';
mygrid.memoryLimit=2;
mygrid.ioSize=[nx*ny 1];
mygrid.facesSize=[nfx' nfy'];
mygrid.facesExpand=[nx 3*nx];
mygrid.missVal=0;

fldstr2d={'XC','YC','XG','YG','RAC','Depth','DXG','DYG','DXC','DYC'};
fldstr3d={'hFacC','hFacW','hFacS','mskC','mskS','mskW'};
fldstr3dp={'hFacC','hFacW','hFacS','maskCtrlC','maskCtrlS','maskCtrlW'};
fldstr1d={'RC','RF','DRC','DRF'};

if(read_netcdf);
  itile=1:29;	%read full aste domain, all 29 tilese
  load([dirRoot '/ASTE_mlibrary/tiles_params.mat']);	%load structure tiles, change absolute path as needed

  grid_nc=nctiles_grid(fldstr1d,dirGrid,itile,tiles,1);	%last flag "1" yields [270 1350] size for all fields
  for i=1:size(fldstr1d,2);mygrid.(fldstr1d{i})=squeeze(grid_nc.(fldstr1d{i}));end;

  grid_nc=nctiles_grid({fldstr3d{1:3}},dirGrid,itile,tiles,1);	%last flag "1" yields [270 1350] size for all fields
  for i=1:3;mygrid.(fldstr3d{i})=squeeze(grid_nc.(fldstr3d{i}));end;
  tmp=mygrid.hFacC;tmp(find(tmp>0))=1;mygrid.mskC=tmp;
  tmp=mygrid.hFacW;tmp(find(tmp>0))=1;mygrid.mskW=tmp;
  tmp=mygrid.hFacS;tmp(find(tmp>0))=1;mygrid.mskS=tmp;

  grid_nc=nctiles_grid(fldstr2d,dirGrid,itile,tiles,1);	%last flag "1" yields [270 1350] size for all fields
  for i=1:size(fldstr2d,2);mygrid.(fldstr2d{i})=squeeze(grid_nc.(fldstr2d{i}));end;

else;
  for i=1:size(fldstr1d,2);tmp=squeeze(rdmds([dirGrid fldstr1d{i}]));mygrid.(fldstr1d{i})=tmp;end;
  for i=1:size(fldstr3d,2);tmp=rdmds([dirGrid fldstr3dp{i}]);
  	tmp=reshape(tmp,nx,ny,nz);mygrid.(fldstr3d{i})=tmp;end;
  for i=1:size(fldstr2d,2);tmp=rdmds([dirGrid fldstr2d{i}]);
          tmp=reshape(tmp,nx,ny);mygrid.(fldstr2d{i})=tmp;end;
end;

mygrid.mskC(find(mygrid.mskC==0))=nan;
							%/*}}}*/ [grid]
clear areaW areaS Vol
for k=1:nz;  areaW(:,:,k)=mygrid.DYG.*mygrid.DRF(k);
             areaS(:,:,k)=mygrid.DXG.*mygrid.DRF(k); Vol(:,:,k)=mygrid.RAC.*mygrid.DRF(k);end;

%mask out obcs:
RAC=mygrid.RAC; if(test3d);  RAC3=repmat(RAC,[1 1 nz]); end;
hfC=mygrid.hFacC;
DD=mygrid.Depth;
dxg=mygrid.DXG;dyg=mygrid.DYG;
dxg3d=repmat(dxg,[1 1 nz]);
dyg3d=repmat(dyg,[1 1 nz]);

hfC(find(hfC==0))=nan;
hfC1=hfC(:,:,1);hfC1(find(hfC1==0))=nan;
RACg=RAC.*hfC1;
hfC1p=get_aste_tracer(hfC1,nfx,nfy);
%from data.obcs:
%% OB_Jsouth =   270*24, 900*0,
%% OB_Ieast  =   260*0, 2*96, 458*0, 270*125, 270*427,
  hfC1p(:,845)=nan;hfC1p(:,24)=nan;hfC1p(366,261:262)=nan;      %Last one: Gibraltar Strait
  hfC1p=aste_tracer2compact(hfC1p,nfx,nfy);
RACgp=RAC.*hfC1p;

istart=1;iend=1;%length(flist)-1;%need -1 to get tendency
if(istart~=1);fprintf('WARNING: starting at: %i\n',istart);end;
if(iend~=length(flist)-1);fprintf('WARNING: ending not at end of list, but at: %i\n',iend);end;

%===== start time loop ===============
for it=istart:iend;
  if1=it;
  if2=if1+1;
  if(read_netcdf);
    t1=flist(if1);
    t2=flist(if2);
  else;
    t1=str2num(flist(if1).name(idot));
    t2=str2num(flist(if2).name(idot));%t1=2232;t2=4248;
  end;
  dt=(t2-t1)*deltaTime;
  t1str=sprintf('%10.10i',t1);t2str=sprintf('%10.10i',t2);

%===================== BLOCK 1: Tendency ======================================================
%ocean
%THETA SALT/*{{{*/
%'UVELMASS' 'VVELMASS' 'ADVx_TH ' 'ADVy_TH ' 'DFxE_TH ' 'DFyE_TH '
%'ADVx_SLT' 'ADVy_SLT' 'DFxE_SLT' 'DFyE_SLT' 'oceSPtnd' 'AB_gT   ' 'AB_gS   '

%flist=dir([dirIn 'budg3d_snap_set2.*.data']);%THETADR,SALTDR
%a=read_slice([dirIn flist(if1).name],nx,ny,1:nz,'real*8');
%b=read_slice([dirIn flist(if2).name],nx,ny,1:nz,'real*8');
%THETA=(b-a)/dt;					%3d, m.degC/s/*}}}*/

  AB_gT=0;%read_slice([dirIn 'budg3d_kpptend_set1.' t2str '.data'],nx,ny,[1:nz]+5*nz,'real*8');

  if(test3d);varName={'THETADR'};n3=nz;else;varName={'THETA'};  n3=1;end;
  if(read_netcdf);
    if(test3d);fileName='nctiles_monthly_snapshots';else;fileName='nctiles_monthly_snapshots_integrated';end;
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if1:if2],tiles,1);
      eval([strtrim(varName{ifld}) '=squeeze(fldOut(:,:,:,2)-fldOut(:,:,:,1))./dt;']);
    end;
  else;
    if(test3d);fileName='budg3d_snap_set2';else;fileName='budg2d_snap_set2';end;
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,[t1,t2]);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=(fldOut.' varName{ifld} '(:,:,:,2)-fldOut.' varName{ifld} '(:,:,:,1))./dt;']);
    end;
  end;
  clear fldOut;   %clear memory
  if exist('THETADR','var');  THETA = THETADR; clear THETADR; end;

  tmptend=myparms.rcp.*(THETA-AB_gT).*mk3D_mod(RAC,THETA);%kg/m3 J/kg/degC degC/s m3 = [J/s]=Watts
  if(test3d); budgO.fluxes.tend=tmptend; end;		%3d Watts
  budgO.tend=nansum(tmptend,3);				%2d Watts

%ice
% 'ETAN    ' 'SIheff  ' 'SIhsnow ' 'SIarea  ' 'sIceLoad' 'PHIBOT  ' %/*{{{*/
%flist=dir([dirIn 'budg2d_snap_set1.*.data']);
%a=readbin([dirIn flist(if1).name],[nx ny 6],1,'real*8');
%b=readbin([dirIn flist(if2).name],[nx ny 6],1,'real*8');
%c=(b-a)./dt;
%SIheff=c(:,:,2);
%SIhsnow=c(:,:,3);/*}}}*/

  varName={'SIheff','SIhsnow'};n3=1;
  if(read_netcdf);
    fileName='nctiles_monthly_snapshots';
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if1:if2],tiles,1);
      eval([strtrim(varName{ifld}) '=squeeze(fldOut(:,:,:,2)-fldOut(:,:,:,1))./dt;']);
    end;
  else;
    fileName='budg2d_snap_set1';
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,[t1,t2]);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=(fldOut.' varName{ifld} '(:,:,:,2)-fldOut.' varName{ifld} '(:,:,:,1))./dt;']);
    end;
  end;
  clear fldOut;   %clear memory

  budgI.tend=-myparms.flami*(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn).*RAC;%J/kg.kg/m3.m/s m2 = J/s = W

%ocean and ice
  budgOI.tend=budgO.tend+budgI.tend;			%2d, %Watts

%===================== BLOCK 2: Vertical convergence ==============================================
%'budg2d_zflux_set1.' t2str '.meta' %/*{{{*/
% 'oceFWflx' 'SIatmFW ' 'TFLUX   ' 'SItflux ' 'SFLUX   ' 'oceQsw  ' 'oceSPflx'
%budg2d_zflux_set2.' t2str '.meta
% 'SRELAX  ' 'TRELAX  ' 'WTHMASS ' 'WSLTMASS' 'oceSflux' 'oceQnet '
% 'SIatmQnt' 'SIaaflux' 'SIsnPrcp' 'SIacSubl'
%/*}}}*/

%this is heat fluxes related to a saltplume scheme that has not been shown to close budget
%so it has been turn off. for now, it should always be zeros.  Read further way below.
  varName={'SPforcT1'};	
  SPforcT1=zeros(nx,ny,1);

  varName={'TFLUX','oceQsw','SItflux'};n3=1;%,'SPforcT1'
  if(read_netcdf);
    fileName='nctiles_monthly';
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if2],tiles,1);
      eval([strtrim(varName{ifld}) '=fldOut(:,:,:,1);']);
    end;
  else;
    fileName='budg2d_zflux_set1';
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,t2);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=fldOut.' varName{ifld} '(:,:,:,1);']);
    end;
  end;
  clear fldOut;

%SIaaflux= , W/m2, conservative ocn<->seaice adv. heat flux adjust, lots of things added to Qnet, in seaice_growth.F/*{{{*/
%%if ~myparms.useRFWF|~myparms.useNLFS;
  %SIaaflux=read_slice([dirIn 'budg2d_zflux_set2.' t2str '.data'],nx,ny,8,'real*8');	%read up there already
%%else;
%%  SIaaflux=0;
%%end;%/*}}}*/

%for ASTE, we do not use any relaxation, so will not distribute all zeros files
%also, WTHMASS is a special case where we write 3d fields to only 2d surface.
  TRELAX=zeros(nx,ny);
  varName={'oceQnet','WTHMASS','SIaaflux'};n3=1;
  if(read_netcdf);
    fileName='nctiles_monthly';
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if2],tiles,1);
      eval([strtrim(varName{ifld}) '=fldOut(:,:,:,1);']);
    end;
  else;
    fileName='budg2d_zflux_set2';
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,t2);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=fldOut.' varName{ifld} '(:,:,:,1);']);
    end;
  end;
  clear fldOut;

%WTHMASS is output as a single layer without k-index, so assume it's the top layer (not depth integr)
%WTHMASS=wVel(i,j,k,bi,bj)*0.5*(theta(i,j,k,bi,bj)+theta(k,j,km1,bi,bj)) <-- between km1 and k, so TOP surface
  if ~myparms.useNLFS;%|myparms.rStar==0;
    %WTHMASS2d=read_slice([dirIn 'budg2d_zflux_set2.' t2str '.data'],nx,ny,3,'real*8');%%C.m/s , read way above
  else;
    WTHMASS=0.*WTHMASS;
  end;

  geothFlux=0;

  if myparms.SPHF
%SPforcT1: heat flux at surface due to SALT_PLUME_VOL (note SPforcS1 = oceSPflx)  %/*{{{*/
%  SPforcT1=read_slice([dirIn 'budg2d_zflux_set1.' t2str '.data'],nx,ny,8,'real*8');%/*}}}*/
%For ASTE R1, we are never in this loop, so no data are written for netcdf/binary
    varName={'oceEPtnd'};if(test3d);n3=nz;else;n3=1;end;
    if(read_netcdf);
      oceEPtnd=0.*ones(nx,ny,nz);
    else;
%      oceEPtnd=read_slice([dirIn 'budg3d_kpptend_set1.' t2str '.data'],nx,ny,[1:nz]+7*nz,'real*8');
      if(test3d);fileName='budg3d_kpptend_set1';else;fileName='budg2d_hflux_set2';end;
      fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,t2);oceEPtnd=fldOut.oceEPtnd;clear fldOut;
    end;
    %SPforcT1: read way above
  else;
    SPforcT1=0.*ones(nx,ny);
    oceEPtnd=0.*ones(nx,ny,nz);
  end;

%Vertical Advective and Diffusive (ex and im) Flux of heat, degC.m^3/s
%skip WTHMASS, only need 1st layer, to be read later; also, for integrated, they are saved 
%in budg2d_kpptend_set1 to reduce number of output files, even though they're not tend
  varName={'ADVr_TH','DFrE_TH','DFrI_TH'};if(test3d);n3=nz;else;n3=1;end;
  if(read_netcdf);
    if(test3d);fileName='nctiles_monthly';else;fileName='nctiles_monthly_integrated';end;
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if2],tiles,1);
      eval([strtrim(varName{ifld}) '=fldOut(:,:,:,1);']);
    end;
  else;
    if(test3d);fileName='budg3d_zflux_set2';else;fileName='budg2d_kpptend_set1';end;
    %WTHMASS=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+1*nz,'real*8');%3d, degC.m/s, 1st layer=WTHMASS2d/*{{{*/
    %ADVr_TH=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+3*nz,'real*8');%3d, degC.m3/s
    %DFrE_TH=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+4*nz,'real*8');%3d, degC.m3/s
    %DFrI_TH=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+5*nz,'real*8');%3d, degC.m3/s/*}}}*/
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,t2);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=fldOut.' varName{ifld} '(:,:,:,1);']);
    end;
  end;
  clear fldOut

%KPP non-local Flux of Pot.Temp; %3d or 2d, degC.m3/s
  varName={'KPPg_TH'};if(test3d);n3=nz;else;n3=1;end;
  if(read_netcdf);
    if(test3d);fileName='nctiles_monthly';else;fileName='nctiles_monthly_integrated';end;
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if2],tiles,1);
      eval([strtrim(varName{ifld}) '=fldOut(:,:,:,1);']);
    end;
  else;
    if(test3d);fileName='budg3d_kpptend_set1';else;fileName='budg2d_kpptend_set1';end;
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,[t2]);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=fldOut.' varName{ifld} '(:,:,:,1);']);
    end;
  end;
  clear fldOut

  budgO.zconv=TFLUX+geothFlux+SPforcT1;				%2d
  budgI.zconv=-(SItflux+TFLUX-TRELAX+SPforcT1);			%2d

%in linear surface we omit :
  if ~myparms.useNLFS;%|myparms.rStar==0; 
    budgO.zconv=budgO.zconv-myparms.rcp*WTHMASS(:,:,kBudget); 	%2d
  end;
  %in virtual salt flux we omit :
  %if ~myparms.useRFWF|~myparms.useNLFS;  %if useRFWF=0 or useNLFS=0
    budgI.zconv=budgI.zconv+SIaaflux; 
  %end;

  budgO.zdia=budgO.zconv;				%2d

%for deep ocean layer :
  if kBudget>1;
    budgO.zconv=-(ADVr_TH+DFrE_TH+DFrI_TH)./RAC*myparms.rcp;%2d (only for layer 11)
    if(myparms.SPHF);budgO.zconv=budgO.zconv+oceEPtnd(:,:,kBudget);end;
    budgO.zdia=-(DFrE_TH+DFrI_TH)./RAC*myparms.rcp;	%2d
    dd=mygrid.RF(kBudget); 
    msk=mygrid.mskC(:,:,kBudget);
    swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
    if dd<-200; swfrac=0; end;
    budgO.zconv=budgO.zconv+swfrac*oceQsw+geothFlux;%.*msk;%2d
    budgO.zdia =budgO.zdia +swfrac*oceQsw+geothFlux;%.*msk;	%2d, w/o advective
  end;
%
%notes: - geothFlux remains to be accounted for in 3D case
%       - diaWtop, diaWbot remain to be implemented in 3D case
  if(test3d);
    nr=size(ADVr_TH,3);
    trWtop=-(ADVr_TH+DFrE_TH+DFrI_TH+KPPg_TH)*myparms.rcp;	%3d, Watts (degC.m3/s kg/m3 J/kg/degC)
    diaWtop=-(DFrE_TH+DFrI_TH+KPPg_TH)*myparms.rcp;		%3d, Watts (degC.m3/s kg/m3 J/kg/degC)
    %
    if(length(mygrid.RF)>nz);%in raw binary
      dd=mygrid.RF(1:end-1);				%50, m
    else;		     %in netcdf
      dd=mygrid.RF(1:end);				%50, m
    end;
    swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);		%unitless
    swfrac(dd<-200)=0;
    swtop=mk3D_mod(swfrac,trWtop).*mk3D_mod(RAC.*oceQsw,trWtop);	%frac sw penetr top, 3d, Watts
    mskC=mygrid.mskC;
    swtop(isnan(mskC))=0;
    trWtop=trWtop+swtop;					%3d,adv+diff+penetr flx, Watts
    diaWtop=diaWtop+swtop;				%3d,    diff+penetr flx, Watts
    %
    if(myparms.SPHF);
      tmp1=(mk3D_mod(SPforcT1,oceEPtnd)-cumsum(oceEPtnd,3)).*RAC3;
      trWtop(:,:,2:nr)=trWtop(:,:,2:nr)+tmp1(:,:,1:nr-1);
      diaWtop(:,:,2:nr)=diaWtop(:,:,2:nr)+tmp1(:,:,1:nr-1);
    end;
    trWtop(:,:,1)=budgO.zconv.*RAC;			%no minus, toplayer = surf flx or layer11
    trWbot(:,:,1:nr-1)=trWtop(:,:,2:length(mygrid.RC));	%3d
    trWbot(:,:,length(mygrid.RC))=0;			%need to take care of geothFlux somewhere
    diaWtop(:,:,1)=budgO.zdia.*RAC;			%no minus, toplayer = surf flx or layer11
    diaWbot(:,:,1:length(mygrid.RC)-1)=diaWtop(:,:,2:length(mygrid.RC));		%3d
    diaWbot(:,:,length(mygrid.RC))=0;			%need to take care of geothFlux somewhere
    %
    budgO.fluxes.trWtop=trWtop;%Watt			%3d
    budgO.fluxes.trWbot=trWbot;%Watt			%3d
    budgO.fluxes.diaWtop=diaWtop;%Watt			%3d
    budgO.fluxes.diaWbot=diaWbot;%Watt			%3d
  else;
    budgO.fluxes.trWtop=-RAC.*(budgO.zconv-geothFlux);	%note [-] sign, not sure why geoth is to the top
    budgO.fluxes.trWbot=RAC.*geothFlux;%Watt		%2d
    budgO.fluxes.diaWtop=-RAC.*(budgO.zdia-geothFlux);	%2d
    budgO.fluxes.diaWbot=RAC.*geothFlux;%Watt		%2d
  end;
  budgI.fluxes.trWtop=-RAC.*(budgI.zconv+budgO.zconv);	%2d
  budgI.fluxes.trWbot=-RAC.*budgO.zconv;%Watt		%2d
%
  budgO.fluxes.zconv=budgO.fluxes.trWtop-budgO.fluxes.trWbot;%3d, Watts	TOP minus BOTTOM
  budgO.zconv=mk3D_mod(RAC,budgO.zconv).*budgO.zconv;	%Watt,%2d, why mk3D?
  budgI.zconv=RAC.*budgI.zconv;%Watt			%2d
  budgOI.zconv=budgO.zconv+budgI.zconv;		%2d

%sign convention: budgO.fluxes.trWtop + buddgI.fluxes.trWbot = identically zero
%so sign means going "INTO" the box, e.g., top = going down, bot = going up

%===================== BLOCK 3: Horizontal convergence ======================================
% budg2d_hflux_set1.0000000336.meta/*{{{*/
% 'ADVxHEFF' 'ADVyHEFF' 'DFxEHEFF' 'DFyEHEFF' 'ADVxSNOW' 'ADVySNOW' 'DFxESNOW' 'DFyESNOW'
% budg2d_hflux_set2.0000000336.meta
% 'UVELMASS' 'VVELMASS' 'ADVx_TH ' 'ADVy_TH ' 'DFxE_TH ' 'DFyE_TH '
% 'ADVx_SLT' 'ADVy_SLT' 'DFxE_SLT' 'DFyE_SLT' 'oceSPtnd' 'AB_gT   ' 'AB_gS   '
%/*}}}*/
  varName={'ADVx_TH','ADVy_TH','DFxE_TH','DFyE_TH'};if(test3d);n3=nz;else;n3=1;end;
  if(read_netcdf);
    if(test3d);fileName='nctiles_monthly';else;fileName='nctiles_monthly_integrated';end;
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if2],tiles,1);
      eval([strtrim(varName{ifld}) '=fldOut(:,:,:,1);']);
    end;
  else;
    if(test3d);fileName='budg3d_hflux_set2';else;fileName='budg2d_hflux_set2';end;
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,[t2]);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=fldOut.' varName{ifld} '(:,:,:,1);']);
    end;
  end;
  clear fldOut

  tmpUo=myparms.rcp*(ADVx_TH+DFxE_TH);		%3d, kg/m3.J/kg/degC.degC.m3/s = [W]
  tmpVo=myparms.rcp*(ADVy_TH+DFyE_TH);		%3d, kg/m3.J/kg/degC.degC.m3/s = [W]

  if(test3d)
      tmp=calc_UV_conv_mod(tmpUo,tmpVo,nfx,nfy);
      budgO.fluxes.hconv=cat(2,tmp{1},tmp{3},reshape(tmp{4},nfy(4),nfx(4),nz),reshape(tmp{5},nfy(5),nfx(5),nz));
  end;

  tmp=calc_UV_conv_mod(nansum(tmpUo,3),nansum(tmpVo,3),nfx,nfy);
  budgO.hconv=cat(2,tmp{1},tmp{3},reshape(tmp{4},nfy(4),nfx(4),1),reshape(tmp{5},nfy(5),nfx(5),1));

  if(test3d);
    tmpUo=myparms.rcp*(ADVx_TH+DFxE_TH);		%3d, kg/m3.J/kg/degC.degC.m3/s = [W]
    tmpVo=myparms.rcp*(ADVy_TH+DFyE_TH);		%3d, kg/m3.J/kg/degC.degC.m3/s = [W]
  end;

%ice, note for ASTER1, all the diff terms DF[x,y]E{HEFF,SNOW} are zeros
  varName={'DFxEHEFF','DFyEHEFF','DFxESNOW','DFyESNOW'};
  for ifld=1:size(varName,2);eval([varName{ifld} '=zeros(nx,ny);']);end;

  varName={'ADVxHEFF','ADVyHEFF','ADVxSNOW','ADVySNOW'};n3=1;
  if(read_netcdf);
    fileName='nctiles_monthly';
    for ifld=1:size(varName,2);
      [fldOut time_step]=nctiles2aste_v2(varName{ifld},[dirIn '/' fileName '/'],itile,1:n3,[if2],tiles,1);
      eval([strtrim(varName{ifld}) '=fldOut(:,:,:,1);']);
    end;
  else;
    fileName='budg2d_hflux_set1';
    fldOut=rdmds_list_fld(fileName,varName,dirIn,nx,ny,n3,[t2]);
    for ifld=1:length(varName);
      eval([strtrim(varName{ifld}) '=fldOut.' varName{ifld} '(:,:,:,1);']);
    end;
  end;
  clear fldOut; %clear memory

%               J/kg        kg/m^3 m.m^2/s         kg/m^3   m.m^2/s = [W]
  tmpUi=-myparms.flami*(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW...
                       +myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
  tmpVi=-myparms.flami*(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW...
                       +myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);

%  budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;                %W

  temp=calc_UV_conv_mod(tmpUi,tmpVi,nfx,nfy);
  budgI.hconv=cat(2,temp{1},temp{3},reshape(temp{4},nfy(4),nfx(4),1),reshape(temp{5},nfy(5),nfx(5),1));

  budgOI.hconv=budgO.hconv+budgI.hconv;                                  %2d, [W]

  if(test3d);
  %               J/kg        kg/m^3 m.m^2/s         kg/m^3   m.m^2/s = [W]
    tmpUi=-myparms.flami*(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW...
                         +myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
    tmpVi=-myparms.flami*(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW...
                         +myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
  end;

%complete list of output same way as in gcmfaces:
  %here we output tendencies and fluxes in W
  %budgHo=budgO; budgHi=budgI; budgHoi=budgOI;
  budgO.specs.units='W';%ocean only
  budgI.specs.units='W';%ice only

%==============
%/*{{{*/
%print some diagnostics to screen
  if(it==1|mod(it,24)==0);
    a=budgOI.tend.*hfC1p;
    b=budgOI.hconv.*hfC1p;
    c=budgOI.zconv.*hfC1p;
    budgetOI(1,1)=nansum(a(:))./nansum(RACgp(:));
    budgetOI(2,1)=nansum(b(:))./nansum(RACgp(:));
    budgetOI(3,1)=nansum(c(:))./nansum(RACgp(:));
    fprintf('%i %6.4e ',[it budgetOI(1)-budgetOI(2)-budgetOI(3)]);
    
    a=budgI.tend.*hfC1p;
    b=budgI.hconv.*hfC1p;
    c=budgI.zconv.*hfC1p;
    budgetI(1,1)=nansum(a(:))./nansum(RACgp(:));
    budgetI(2,1)=nansum(b(:))./nansum(RACgp(:));
    budgetI(3,1)=nansum(c(:))./nansum(RACgp(:));
    fprintf('%6.4e ',budgetI(1)-budgetI(2)-budgetI(3));
    
    a=budgO.tend.*hfC1p;
    b=budgO.hconv.*hfC1p;
    c=budgO.zconv.*hfC1p;
    budgetO(1,1)=nansum(a(:))./nansum(RACgp(:));
    budgetO(2,1)=nansum(b(:))./nansum(RACgp(:));
    budgetO(3,1)=nansum(c(:))./nansum(RACgp(:));
    fprintf('%6.4e ',budgetO(1)-budgetO(2)-budgetO(3));
  end;

  %if(test3d)
    %a=budgO.fluxes.tend.*repmat(hfC1p,[1 1 nz]);
    %b=budgO.fluxes.hconv.*repmat(hfC1p,[1 1 nz]);
    %c=budgO.fluxes.zconv.*repmat(hfC1p,[1 1 nz]);
    %a(isnan(a))=0;b(isnan(b))=0;c(isnan(c))=0;
  %end;%test3d
%/*}}}*/

%============================================
%now sum up globally
  if(save_budg_scalar==1);		%/*{{{*/ [save_budg_scalar]

%--------------------
%create a list:
    listDiags=['glo_heat_ocn glo_heat_tot glo_heat_ice'];	%/*{{{*/ [2]
%reformat listDiags to cell array:
    jj=strfind(listDiags,' '); jj=[0 jj length(listDiags)+1];
    for ii=1:length(jj)-1;
      tmp1=listDiags(jj(ii)+1:jj(ii+1)-1);
      if ii==1; listDiags2={tmp1}; else; listDiags2{ii}=tmp1; end;
    end;
    listDiags=listDiags2; clear listDiags2;

    if(save_budg_lev==1);
      listDiagsLev{1}='glo_heat_tot_lev';
    end;							%/*}}}*/ [2]
%--------------------

%now do global sum 
%after this sum, unit is: [W/m2]				/*{{{*/ [2]
    tmp1=calc_mskmean_T_mod(budgOI,hfC1p,RACgp,'extensive');glo_heat_tot=[tmp1.tend;tmp1.hconv;tmp1.zconv];
    tmp1=calc_mskmean_T_mod(budgO,hfC1p,RACgp,'extensive'); glo_heat_ocn=[tmp1.tend;tmp1.hconv;tmp1.zconv];
    tmp1=calc_mskmean_T_mod(budgI,hfC1p,RACgp,'extensive'); glo_heat_ice=[tmp1.tend;tmp1.hconv;tmp1.zconv];

%diags_select.m
  %package results
    onediag=[];
    onediag.listTimes=myparms.yearFirst(1)+t1*myparms.timeStep/86400/365.25;
    onediag.listSteps=t1;
    for jj=1:length(listDiags);
        eval(['onediag.' listDiags{jj} '=' listDiags{jj} ';']);
    end;
								%/*}}}*/ [2]
%============================================
    if(save_budg_lev==1);							%/*{{{*/ [2]
      budg3d=budgO.fluxes;
      budg3d.tend(:,:,1)=budg3d.tend(:,:,1)+budgI.tend;
      budg3d.hconv(:,:,1)=budg3d.hconv(:,:,1)+budgI.hconv;
      budg3d.zconv(:,:,1)=budg3d.zconv(:,:,1)+budgI.zconv;
      for k=1:size(kz,1);
        iz=kz(k,1):kz(k,2);
        budglev.tend =squeeze(nansum(budg3d.tend(:,:,iz).*hfC(:,:,iz),3)).*hfC1p;
        budglev.hconv=squeeze(nansum(budg3d.hconv(:,:,iz).*hfC(:,:,iz),3)).*hfC1p;
        budglev.zconv=squeeze(nansum(budg3d.zconv(:,:,iz).*hfC(:,:,iz),3)).*hfC1p;
        tmp1=calc_mskmean_T_mod(budglev,hfC1p,RACgp,'extensive'); glo_heat_tot_lev(k,:)=[tmp1.tend;tmp1.hconv;tmp1.zconv];
      end;

      for jj=1:length(listDiagsLev);
        eval(['onediag.' listDiagsLev{jj} '=' listDiagsLev{jj} ';']);
      end;

    end;									%/*}}}*/ [2]

    fileMat=['diags_set_D_' sprintf('%10.10i',t2) '_' lower(strbudg) '.mat'];%bug-fix: changed t1-->t2
    diags_store_to_mat(dirMat,fileMat,onediag);
  end;%save_budg_scalar			%/*}}}*/ [save_budg_scalar]
%============================================

%plot some figures to show if budgets are closed at every single cell:
  if(plot_fig==1);
    if(it==1|it==length(flist)-1);
      msk1=get_aste_tracer(hfC1p,nfx,nfy);msk1(1:nx,nfy(1)+1:nfy(1)+nx+nfx(4))=nan;

      if(test3d) %/*{{{*/
        figure(1);clf;
        aa=zeros(50,4);
        n=1;pt(n,:)=floor([17 226]./fac);ix=pt(n,1);iy=pt(n,2);k=5;
        bb =[budgO.fluxes.tend(ix,iy,k) budgO.fluxes.zconv(ix,iy,k) budgO.fluxes.hconv(ix,iy,k)];
        bbp=[budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)];
        k=1:nz;aa(:,n)=squeeze([budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)]);
        subplot(2,2,n);plot(aa(:,n),-[1:50],'.-');grid;xlabel(['net budg' strbudg(1) 'o']);
          title(['[ix,iy]=[' num2str(pt(n,1)) ',' num2str(pt(n,2)) ']; ' ...
                  num2str(max(abs(aa(:,n)))/max(abs(budgO.fluxes.zconv(ix,iy,:))).*100) '%']);
        
        n=2;pt(n,:)=floor([63 92+nfy(1)]./fac);ix=pt(n,1);iy=pt(n,2);k=1;%a very bad point (for Heat) in the ARctic
        bb =[budgO.fluxes.tend(ix,iy,k) budgO.fluxes.zconv(ix,iy,k) budgO.fluxes.hconv(ix,iy,k)];
        bbp=[budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)];
        k=1:nz;aa(:,n)=squeeze([budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)]);
        subplot(2,2,n);plot(aa(:,n),-[1:50],'.-');grid;xlabel(['net budg' strbudg(1) 'o']);
          title(['[ix,iy]=[' num2str(pt(n,1)) ',' num2str(pt(n,2)) ']; ' ...
                  num2str(max(abs(aa(:,n)))/max(abs(budgO.fluxes.zconv(ix,iy,:))).*100) '%']);

        n=3;pt(n,:)=floor([117 37+nfy(1)]./fac);ix=pt(n,1);iy=pt(n,2);k=1;
        bb =[budgO.fluxes.tend(ix,iy,k) budgO.fluxes.zconv(ix,iy,k) budgO.fluxes.hconv(ix,iy,k)];
        bbp=[budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)];
        k=1:nz;aa(:,n)=squeeze([budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)]);
        subplot(2,2,n);plot(aa(:,n),-[1:50],'.-');grid;xlabel(['net budg' strbudg(1) 'o']);
          title(['[ix,iy]=[' num2str(pt(n,1)) ',' num2str(pt(n,2)) ']; ' ...
                  num2str(max(abs(aa(:,n)))/max(abs(budgO.fluxes.zconv(ix,iy,:))).*100) '%']);

        n=4;pt(n,:)=floor([88 107+nfy(1)]./fac);ix=pt(n,1);iy=pt(n,2);k=1;
        bb =[budgO.fluxes.tend(ix,iy,k) budgO.fluxes.zconv(ix,iy,k) budgO.fluxes.hconv(ix,iy,k)];
        bbp=[budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)];
        k=1:nz;aa(:,n)=squeeze([budgO.fluxes.tend(ix,iy,k)-budgO.fluxes.zconv(ix,iy,k)-budgO.fluxes.hconv(ix,iy,k)]);
        subplot(224);plot(aa(:,n),-[1:50],'.-');grid;xlabel(['net budg' strbudg(1) 'o']);
          title(['[ix,iy]=[' num2str(pt(n,1)) ',' num2str(pt(n,2)) ']; ' ...
                  num2str(max(abs(aa(:,n)))/max(abs(budgO.fluxes.zconv(ix,iy,:))).*100) '%']);

        figure(1);set(gcf,'paperunit','inches','paperposition',[0 0 9.5 8.5]);
        fpr=[dirOut strbudg '_budget1' '_' sprintf('%10.10i',t1) '.png'];print(fpr,'-dpng');
      end;
%/*}}}*/
      if(test3d);%/*{{{*/
        a=get_aste_tracer(budgO.fluxes.tend.*repmat(hfC1p,[1 1 nz]),nfx,nfy);
        b=get_aste_tracer(budgO.fluxes.hconv.*repmat(hfC1p,[1 1 nz]),nfx,nfy);
        c=get_aste_tracer(budgO.fluxes.zconv.*repmat(hfC1p,[1 1 nz]),nfx,nfy);
        d=get_aste_tracer(swtop.*repmat(hfC1p,[1 1 nz]),nfx,nfy);
        a(isnan(a))=0;b(isnan(b))=0;c(isnan(c))=0;d(isnan(d))=0;
        msk=get_aste_tracer(hfC1p,nfx,nfy);msk(1:nx,nfy(1)+1:nfy(1)+nx+nfx(4))=nan;
        klev=[1,2,41];
        figure(2);clf;colormap(seismic(21));
        for k=1:length(klev);
          k1=klev(k);
          str=['; k=' num2str(k1)];
          subplot(3,4,(k-1)*4+1);mypcolor(a(:,:,k1)'.*(1-isnan(msk1))');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
                 shadeland(1:2*nx,1:nfy(1)+nx+nfx(4),isnan(msk1'),[.7 .7 .7]);title([strbudg ' tend' str]);
          subplot(3,4,(k-1)*4+2);mypcolor(b(:,:,k1)'.*(1-isnan(msk1))');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
                 shadeland(1:2*nx,1:nfy(1)+nx+nfx(4),isnan(msk1'),[.7 .7 .7]);title([strbudg ' hconv']);
          subplot(3,4,(k-1)*4+3);mypcolor(c(:,:,k1)'.*(1-isnan(msk1))');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
                 shadeland(1:2*nx,1:nfy(1)+nx+nfx(4),isnan(msk1'),[.7 .7 .7]);title([strbudg ' zconv']);
          subplot(3,4,(k-1)*4+4);mypcolor((a(:,:,k1)'-b(:,:,k1)'-c(:,:,k1)').*(1-isnan(msk1))');cc1=0.05*max(abs(caxis));caxis([-cc1 cc1]);
                 mythincolorbar;grid;shadeland(1:2*nx,1:nfy(1)+nx+nfx(4),isnan(msk1'),[.7 .7 .7]);title([strbudg ' tend-hconv-zconv']);
        end;
        figure(2);set(gcf,'paperunit','inches','paperposition',[0 0 14 9]);
        fpr=[dirOut strbudg '_budget2' '_' sprintf('%10.10i',t1) '.png'];print(fpr,'-dpng');
      end;%/*}}}*/

      if(test3d) %/*{{{*/
        a=budgO.fluxes.tend.*repmat(hfC1p,[1 1 nz]);
        b=budgO.fluxes.hconv.*repmat(hfC1p,[1 1 nz]);
        c=budgO.fluxes.zconv.*repmat(hfC1p,[1 1 nz]);
        a(isnan(a))=0;b(isnan(b))=0;c(isnan(c))=0;
        klev=[1,2,41];
        ix=floor(76./fac):nx;iy=nfy(1)+1:nfy(1)+floor(200/fac);
        figure(3);clf;colormap(seismic(21));
        for k=1:length(klev);
          k1=klev(k);
          str=['; k=' num2str(k1)];
          subplot(3,4,(k-1)*4+1);mypcolor(ix,iy,a(ix,iy,k1)');cc1=0.04*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
                 title([strbudg ' tend' str]);shadeland(ix,iy,isnan(hfC1p(ix,iy)'),[.7 .7 .7]);
          subplot(3,4,(k-1)*4+2);mypcolor(ix,iy,b(ix,iy,k1)');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
                 title([strbudg ' hconv']);shadeland(ix,iy,isnan(hfC1p(ix,iy)'),[.7 .7 .7]);
          subplot(3,4,(k-1)*4+3);mypcolor(ix,iy,c(ix,iy,k1)');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
                 title([strbudg ' zconv']);shadeland(ix,iy,isnan(hfC1p(ix,iy)'),[.7 .7 .7]);
          subplot(3,4,(k-1)*4+4);mypcolor(ix,iy,a(ix,iy,k1)'-b(ix,iy,k1)'-c(ix,iy,k1)');cc1=0.05*max(abs(caxis));caxis([-cc1 cc1]);
                 mythincolorbar;grid;title([strbudg ' tend-hconv-zconv']);shadeland(ix,iy,isnan(hfC1p(ix,iy)'),[.7 .7 .7]);
        end;
        figure(3);set(gcf,'paperunit','inches','paperposition',[0 0 14 9]);
        fpr=[dirOut strbudg '_budget3' '_' sprintf('%10.10i',t1) '.png'];print(fpr,'-dpng');
      end;	%test3d
%/*}}}*/
      a=get_aste_tracer(budgO.tend.*hfC1p ,nfx,nfy); d=get_aste_tracer(budgI.tend.*hfC1p,nfx,nfy);%/*{{{*/
      b=get_aste_tracer(budgO.hconv.*hfC1p,nfx,nfy); f=get_aste_tracer(budgI.hconv.*hfC1p,nfx,nfy);
      c=get_aste_tracer(budgO.zconv.*hfC1p,nfx,nfy); g=get_aste_tracer(budgI.zconv.*hfC1p,nfx,nfy);
      a(isnan(a))=0;b(isnan(b))=0;c(isnan(c))=0;
      d(isnan(d))=0;f(isnan(f))=0;g(isnan(g))=0;
      ix=75:2*nx;iy=1:870;
      k=1;
      str=['; k=' num2str(k)];
      figure(4);clf;colormap(seismic(21));
      subplot(241);mypcolor(ix,iy,a(ix,iy,k)');cc1=0.04*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);title(['ocn ' strbudg ' tend, all z']);
      subplot(242);mypcolor(ix,iy,b(ix,iy,k)');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);title(['ocn ' strbudg ' hconv, all z']);
      subplot(243);mypcolor(ix,iy,c(ix,iy,k)');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);title(['ocn ' strbudg ' zconv, all z']);
      subplot(244);mypcolor(ix,iy,a(ix,iy,k)'-b(ix,iy,k)'-c(ix,iy,k)');cc1=0.05*max(abs(caxis));caxis([-cc1 cc1]);
                   mythincolorbar;grid;title('ocean tend-hconv-zconv, all z');
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);
      subplot(245);mypcolor(ix,iy,d(ix,iy,k)');cc1=0.04*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);title(['ice ' strbudg ' tend']);
      subplot(246);mypcolor(ix,iy,f(ix,iy,k)');cc1=0.1*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid;
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);title(['ice ' strbudg ' hconv']);
      subplot(247);mypcolor(ix,iy,g(ix,iy,k)');cc1=0.04*max(abs(caxis));caxis([-cc1 cc1]);mythincolorbar;grid
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);title(['ice ' strbudg ' zconv']);
      subplot(248);mypcolor(ix,iy,d(ix,iy,k)'-f(ix,iy,k)'-g(ix,iy,k)');cc1=0.05*max(abs(caxis));caxis([-cc1 cc1]);
                   mythincolorbar;grid;title('ice tend-hconv-zconv');
             shadeland(ix,iy,isnan(msk1(ix,iy)'),[.7 .7 .7]);
      figure(4);set(gcf,'paperunit','inches','paperposition',[0 0 14 9]);
      fpr=[dirOut strbudg '_budget4' '_' sprintf('%10.10i',t1) '.png'];print(fpr,'-dpng');
%/*}}}*/
    end;
%keyboard
  end;%plot_fig

%============================================
%saving budgets in 3D, need large storage space!
  if(save_budg_3d);		%/*{{{*/
    tmpbudgO=budgO;
    clear budgO;
    budgO.fluxes=tmpbudgO.fluxes;
    budgO.fluxes=rmfield(budgO.fluxes,'zconv');
    fileMat=['diags_set_D_' sprintf('%10.10i',t2) '_' lower(strbudg) '_3d.mat'];%bug-fix: change t1-->t2
    save([dirMat '/' fileMat],'budgO');
    %save([dirMat '/' fileMat],'budgO','budgI');
  end;				%/*}}}*/
  if(save_budg_2d);		%/*{{{*/
    if(save_budg_3d);budgO=tmpbudgO;end;
    budgO=rmfield(budgO,'fluxes');
    fileMat=['diags_set_D_' sprintf('%10.10i',t2) '_' lower(strbudg) '_2d.mat'];%bug-fix: change t1-->t2
    save([dirMat '/' fileMat],'budgO','budgI');
  end;				%/*}}}*/
%============================================
%conserve outside sea ice area!!!!
%fix: add SIaaflux for when NLFS=2 and RFWF=1

%%load(['diags_set_D_' num2str(t2) '.mat'],...
%load([dirIn '../../mat/diags_set_D/diags_set_D_' num2str(t2) '.mat'],...
%      'glo_heat_ice','glo_heat_ocn','glo_heat_tot');
%glo_heat_tot-budgetOI			%[0.288657986402541  -0.124384150538271   0.612843109593086] x 10^-13
%
end;%ts loop
fprintf('\n');
