%atn 16.jan.2021
%script to mimic eccov4r3 salt 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='Salt';
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.seaice_variable_salinity_flag=0;	%so that we can overwrite myparms.SIsal0
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   '
%note: if full ocean depth, SALT=sum_k(dSALT(:,:,k)/dt*drF(k).*hFacC(:,:,k)) [psu m/s]
%      if(test3d): SALT(k) = SALTDR(k) = dSALT(:,:,k)/dt*drF(k).*hFacC(:,:,k) <-- new diag atn put in
%/*}}}*/
%/*{{{*/

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

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

  if(test3d);varName={'SALTDR'};n3=nz;else;varName={'SALT'};  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('SALTDR','var');  SALT = SALTDR; clear SALTDR; end;

  tmptend=myparms.rhoconst.*(SALT-AB_gS).*mk3D_mod(RAC,SALT);%kg/m3 g/kg.m/s m2 = [g/s]
  if(test3d); budgO.fluxes.tend=tmptend; end;		%3d g/s
  budgO.tend=nansum(tmptend,3);				%2d g/s

%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);
%/*}}}*/
  varName={'ETAN','SIheff'};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
  if(myparms.seaice_variable_salinity_flag==0);%constant salinity case; beginning iter10 we are here.
    budgI.tend=myparms.SIsal0*SIheff*myparms.rhoi.*RAC;%g/kg m/s kg/m3 m2 = g/s
  else;
%this is not correct, will need to figure out what to do with SEAICE_saltFrac and SSS
%but for now, we're not in here, so set to 0.  Will need revisit if end up here.
    budgI.tend=0.*SIheff*myparms.rhoi.*RAC;%g/kg m/s kg/m3 m2 = g/s
  end;

%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'

%2d fluxes:

%total salt flux (match salt-content variations), >0 increases salt, g/m^2/s
%SFLUX = surfaceForcingS(i,j,bi,bj)*rhoConst; (rUnit2mass=rhoConst), surfaceForcingS = surForcS
% if nonlinFreeSurf & useRealFreshWaterFlux: SFLUX=SFLUX+PmEpR.*salt(i,j,1,bi,bj)
%SFLUX   =read_slice([dirIn 'budg2d_zflux_set1.' t2str '.data'],nx,ny,5,'real*8');%2d, g/m2/s
%net surface Salt flux rejected into the ocean during freezing, (+=down), g/m^2/s
%oceSPflx = saltPlumeFlux
%oceSPflx=read_slice([dirIn 'budg2d_zflux_set1.' t2str '.data'],nx,ny,7,'real*8');%2d, g/m2/s
%oceSPflx=0;	%looks like oceSPflx is identically -SFLUX, this makes vertical FW flux zero, wrong!
%SRELAX  =read_slice([dirIn 'budg2d_zflux_set2.' t2str '.data'],nx,ny,1,'real*8');%2d, all zeros
%salt tendency due to salt plume flux >0 increases salinity, g/m^2/s
%oceSPtnd=read_slice([dirIn 'budg3d_hflux_set2.' t2str '.data'],nx,ny,[1:nz]+10*nz,'real*8');%2d,g/m2/s
%/*}}}*/

  varName={'SFLUX','oceSPflx'}; 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_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;

%for ASTE, we do not use any relaxation, so will not distribute all zeros files
%also, WSLTMASS is a special case where we write 3d fields to only 2d surface.
  SRELAX=zeros(nx,ny);
  varName={'oceSflux','WSLTMASS'};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;

  varName={'oceSPtnd'};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;
    %oceSPtnd=read_slice([dirIn 'budg3d_kpptend_set1.' t2str '.data'],nx,ny,[1:nz]+4*nz,'real*8');%3d,g/m2/s
    %oceSPtnd=read_slice([dirIn 'budg2d_hflux_set2.' t2str '.data'],nx,ny,[1:1]+10*1,'real*8');%2d,g/m2/s
    if(test3d);fileName='budg3d_kpptend_set1';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;
%/*{{{*/
%#ifndef SALT_PLUME_VOLUME
%          plumefrac(i,j) = ( plumekb2D(i,j)-plumefrac(i,j) )
%     &                     *maskC(i,j,k,bi,bj)
%          plumetend(I,J) = saltPlumeFlux(i,j,bi,bj)*plumefrac(I,J)
%#else /* SALT_PLUME_VOLUME */
%          plumetend(I,J) = SPforcingS(i,j,k,bi,bj)
%#endif /* SALT_PLUME_VOLUME */
%/*}}}*/
%working approach for real fresh water (?) and virtual salt flux
  if ~myparms.useRFWF|~myparms.useNLFS;	%enters if useRFWF=0 or useNLFS=0
% net surface Salt flux into the ocean (+=down), >0 increases salinity, g/m^2/s
  %oceSflux=read_slice([dirIn 'budg2d_zflux_set2.' t2str '.data'],nx,ny,5,'real*8');%2d, g/m2/s
  %oceSflux: already read above
     error('code missing! need to work out how to close budget in this loop');
  else;
    oceSflux=0.*oceSflux;
  end;

%Vertical Advective and Diffusive (ex and im) Flux of Salinity, psu.m^3/s
%skip WSLTMASS, 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_SLT','DFrE_SLT','DFrI_SLT'};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;
  %ADVr_SLT=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+6*nz,'real*8');%3d, psu.m3/s/*{{{*/
  %DFrE_SLT=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+7*nz,'real*8');%3d, psu.m3/s
  %DFrI_SLT=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+8*nz,'real*8');%3d, psu.m3/s%/*}}}*/
    if(test3d);fileName='budg3d_zflux_set2';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

%KPP non-local Flux of Salt
  varName={'KPPg_SLT'};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	%3d or 2d, psu.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

%in linear surface we omit :
  if ~myparms.useNLFS;
%Vertical Mass-Weight Transp of Salinity, psu.m/s/*{{{*/
%Note: there is NO depth integr for WSLTMASS, and it's read in as 2d field, so must have been layer 1
%      If it's layer 11, Gael added a Depth layer = 11 field, so this is a single layer quantity
  %WSLTMASS2d=read_slice([dirIn 'budg2d_zflux_set2.' t2str '.data'],nx,ny,4,'real*8');	%2d
  %WSLTMASS=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+2*nz,'real*8');%3d, psu.m/s
  %tmp=WSLTMASS2d-WSLTMASS(:,:,1);
  %fprintf('sum(WSLTMASS2d-WSLTMASS(:,:,1)): %10.8e\n',sum(tmp(:)));
  %already read in above, so in this loop we do nothing /*}}}*/
  else;
    WSLTMASS=0.*WSLTMASS;
  end;

  budgO.zconv=SFLUX+oceSPflx;				%2d
  budgI.zconv=-budgO.zconv+SRELAX;			%2d
  %in linear surface we omit :
  if ~myparms.useNLFS;
    budgO.zconv=budgO.zconv-myparms.rhoconst*WSLTMASS(:,:,kBudget); 	%2d
  end;
%in virtual salt flux we omit :
  if ~myparms.useRFWF|~myparms.useNLFS; budgI.zconv=-oceSflux; end;%problematic: what if useRFWF and ~useNLFS?
%atn: attempt1: remove useNLFS filter: no, this is wrong
%if ~myparms.useRFWF; budgI.zconv=-oceSflux; end;
%
  budgO.zdia=budgO.zconv;				%2d

%for deep ocean layer :
  if kBudget>1;
    budgO.zconv=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)./RAC*myparms.rhoconst;%2d (only for layer 11)
    budgO.zconv=budgO.zconv+oceSPtnd(:,:,kBudget);	%2d, why plus here and minus in test3d block?
    budgO.zdia=-(DFrE_SLT+DFrI_SLT)./RAC*myparms.rhoconst;
    budgO.zdia=budgO.zdia+oceSPtnd(:,:,kBudget);%.*msk;
  end;
%
  if(test3d);
    nr=length(mygrid.RC);
    trWtop=-(ADVr_SLT+DFrE_SLT+DFrI_SLT+KPPg_SLT)*myparms.rhoconst;	%3d, g/s
    tmp1=mk3D_mod(oceSPflx,oceSPtnd)-cumsum(oceSPtnd,3);
    tmp1=tmp1.*RAC3;
    trWtop(:,:,2:nr)=trWtop(:,:,2:nr)+tmp1(:,:,1:nr-1);
  % if(~myparms.useNLFS);
  %  trWtop(:,:,1)=(oceSPflx-myparms.rhoconst.*WSLTMASS(:,:,1)).*RAC;	%no minus, toplayer = surf flx or layer11
  % else;
    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
    %
    budgO.fluxes.trWtop=trWtop;%g/s			%3d
    budgO.fluxes.trWbot=trWbot;%g/s			%3d
  else;
    budgO.fluxes.trWtop=-RAC.*budgO.zconv;		%note [-] sign , g/s
    budgO.fluxes.trWbot=RAC.*0;%				%2d, g/s
    budgO.fluxes.diaWtop=-RAC.*budgO.zdia;		%2d
    budgO.fluxes.diaWbot=RAC.*0;%g/s			%2d
  end;
  budgI.fluxes.trWtop=-RAC.*0;				%2d
  budgI.fluxes.trWbot=budgO.fluxes.trWtop(:,:,1);	%2d, g/s
%
  budgO.fluxes.zconv=budgO.fluxes.trWtop-budgO.fluxes.trWbot;%3d, g/s
  budgO.zconv=mk3D_mod(RAC,budgO.zconv).*budgO.zconv;	%g/s,%2d, why mk3D?
  budgI.zconv=RAC.*budgI.zconv;%g/s			%2d
  budgOI.zconv=budgO.zconv+budgI.zconv;		%2d

%===================== 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   '
%ADVx_SLT=read_slice([dirIn 'budg3d_hflux_set2.' t2str '.data'],nx,ny,[1:nz]+6*nz,'real*8');
%ADVy_SLT=read_slice([dirIn 'budg3d_hflux_set2.' t2str '.data'],nx,ny,[1:nz]+7*nz,'real*8');
%DFxE_SLT=read_slice([dirIn 'budg3d_hflux_set2.' t2str '.data'],nx,ny,[1:nz]+8*nz,'real*8');
%DFyE_SLT=read_slice([dirIn 'budg3d_hflux_set2.' t2str '.data'],nx,ny,[1:nz]+9*nz,'real*8');
%/*}}}*/

  varName={'ADVx_SLT','ADVy_SLT','DFxE_SLT','DFyE_SLT'};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.rhoconst*(ADVx_SLT+DFxE_SLT);		%2d or 3d, kg/m3.g/kg.m3/s = [g/s]
  tmpVo=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);		%2d or 3d, kg/m3.g/kg.m3/s = [g/s]

  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.rhoconst*(ADVx_SLT+DFxE_SLT);		%2d or 3d, kg/m3.g/kg.m3/s = [g/s]
    tmpVo=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);		%2d or 3d, kg/m3.g/kg.m3/s = [g/s]
  end;

%ice
%all diffusion terms are zeros for ice and snow in ASTE R1
  DFxEHEFF=zeros(nx,ny);
  DFyEHEFF=zeros(nx,ny);
  varName={'ADVxHEFF','ADVyHEFF'};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

%
  if(myparms.seaice_variable_salinity_flag==0);
    SIsal0=myparms.SIsal0;
  else;
    SIsal0=0;	%not true, would need to calc from SEAICEsaltFrac*SSS, but by accident my prev runs SIsal0=0;
  end;
  tmpUi=SIsal0*myparms.rhoi*(DFxEHEFF+ADVxHEFF);	%g/kg.kg/m3.m.m2/s = g/s
  tmpVi=SIsal0*myparms.rhoi*(DFyEHEFF+ADVyHEFF);
                   
  tmp=calc_UV_conv_mod(tmpUi,tmpVi,nfx,nfy);
  budgI.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);
  tmpUi=SIsal0*myparms.rhoi*(DFxEHEFF+ADVxHEFF);	%g/kg.kg/m3.m.m2/s = g/s
  tmpVi=SIsal0*myparms.rhoi*(DFyEHEFF+ADVyHEFF);
  end;

  budgOI.hconv=budgO.hconv+budgI.hconv;                                  %2d, [g/s]

%complete list of output same way as in gcmfaces:
  %here we output tendencies and fluxes in g/s

  budgO.specs.units='g/s';%ocean only
  budgI.specs.units='g/s';%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;
%/*}}}*/

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

%--------------------
%create a list:
    listDiags=['glo_salt_ocn glo_salt_tot glo_salt_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_salt_tot_lev';
    end;							%/*}}}*/ [2]
%--------------------

%now do global sum 
%after this sum, unit is: [g/s/m2]				/*{{{*/ [2]
    tmp1=calc_mskmean_T_mod(budgOI,hfC1p,RACgp,'extensive');glo_salt_tot=[tmp1.tend;tmp1.hconv;tmp1.zconv];
    tmp1=calc_mskmean_T_mod(budgO,hfC1p,RACgp,'extensive'); glo_salt_ocn=[tmp1.tend;tmp1.hconv;tmp1.zconv];
    tmp1=calc_mskmean_T_mod(budgI,hfC1p,RACgp,'extensive'); glo_salt_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_salt_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);
        
        a(isnan(a))=0;b(isnan(b))=0;c(isnan(c))=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.1*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;				%/*}}}*/
%============================================

end;%ts loop
fprintf('\n');
