%atn 16.jan.2021
%script to mimic eccov4r3 mass 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='Mass';
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]);
drf3d=mk3D_mod(mygrid.DRF,zeros(nx,ny,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
  varName={'ETAN','PHIBOT'}; 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(test3d);
  %3D, with rStar:
    if(myparms.useNLFS<2|myparms.rStar==0);
      tmptend=zeros(nx,ny,nz);%3d
      tmptend(:,:,1)=ETAN*myparms.rhoconst.*RAC; %??%atn 5.sep.2018: this is true only when useRFWF=.true.
  %atn 5.sep.2018 add wrapper for RFWF=0
      if(~myparms.useRFWF);
        tmptend(:,:,1)=zeros(nx,ny);	%for linearFS and useRFWF=.false, the volume is not changing
      end;
  %5.sep.2018
    else;
      if(myparms.useRFWF);
        tmp1=mk3D_mod(mygrid.DRF,hfC).*hfC;					%3d
        tmp2=tmp1./mk3D_mod(DD,tmp1);						%3d
      else;%14.jan.2019: atn: not sure if this else is correct, need to check w/ run when useRFWF=.false. 
        tmp2=drf3d./mk3D_mod(DD,tmp1);					%3d
      end;
      %dETAN/dt*fraction dz/Depth*hFac*rhoconst 				%[m/s kg/m3]
      tmptend=tmp2.*mk3D_mod(ETAN,tmp2)*myparms.rhoconst.*mk3D_mod(RAC,hfC);	%3d, kg/s
    end;
  else;
    tmptend=myparms.rhoconst*ETAN.*RAC;		%2d, [kg/m3 * m/s] = kg/s/m2 * m2
  end;

%mass/second = myparms.rhoconst * sea level / second * RAC
  bp=myparms.rhoconst/9.81*PHIBOT.*RAC;
  if(test3d); budgO.fluxes.tend=tmptend; end;           %3d, still need to work this out
  budgO.tend=nansum(tmptend,3);                        %2d kg/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);ETAN   =c(:,:,1);
%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=(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn).*RAC; %2d, [m/s * kg/m3 *m2]=kg/s

%ocean and ice
  budgOI.tend=budgO.tend+budgI.tend;			%2d, %kg/s

%===================== BLOCK 2: Vertical convergence ==============================================
%'budg2d_zflux_set1.' t2str '.meta' %/*{{{*/
% 'oceFWflx' 'SIatmFW ' 'TFLUX   ' 'SItflux ' 'SFLUX   ' 'oceQsw  ' 'oceSPflx'
%budg3d_zflux_set2.' t2str '.meta
% 'WVELMASS'
%SIatmFW = F(A->I) + F(A->O)
%SIempmr = F(I->O) + F(A->O) = EmPmR
%oceFWflx= - (F(I->O) + F(A->O)) = -SIempmr = -EmPmR
%SIsnPrcp= F(A->I) 
%(SIatmFW - oceFWflx) - SIsnPrcp = F(A->I)+F(A->O)-F(A->O)-F(I->O)-F(A->I) = -F(I->O) = F(O->I)
%/*}}}*/

  varName={'SIsnPrcp'};n3=1;	%kg/m2/s , now precip. (+=dw) over Sea-Ice (area weighted)
  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={'oceFWflx','SIatmFW'};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;
  
  if(test3d)
    varName={'WVELMASS'};n3=nz;
    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;
      %WVELMASS=read_slice([dirIn 'budg3d_zflux_set2.' t2str '.data'],nx,ny,[1:nz]+0*nz,'real*8');%m/s
      fileName='budg3d_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;
  end;
  budgO.zconv=oceFWflx;		%2d, net FW into ocean , +=down, >0 = S decreases , [kg/m2/s]
  budgI.zconv=SIatmFW-oceFWflx;	%2d, net FW minus ocean = on ice

%in virtual salt flux we omit :
  if ~myparms.useRFWF; budgO.zconv=0*budgO.zconv; end	%2d
%for deep ocean with ONLY ONE layer: (need WVELMASS to be 2d)
%note: base on this, because the bottom is zero, so only take top
  if kBudget>1; budgO.zconv=-WVELMASS(:,:,kBudget)*myparms.rhoconst; end; %2d (Gael output only layer 11, not depth integr)
%for deep ocean layer , 3d :
  if(test3d)
    trWtop=-WVELMASS.*myparms.rhoconst;		%3d, note this seems WVELMASS is applied to upper, not lower
    if(myparms.useNLFS<2|myparms.rStar==0);	%atn wrapper
%this is commented out by Gael, why?, atn: add in.  O(10^-3)
%essentially need to remove WVELMASS for the top layer (O(1))
      trWtop(:,:,1)=oceFWflx;			
%atn 5.sep.2018: add another wrapper in attempt to close budget for RFWF=0
      if(~myparms.useRFWF);
        trWtop(:,:,1)=-WVELMASS(:,:,1).*myparms.rhoconst;
      end;
%end 5.sep.2018
   end;						%atn wrapper

    trWbot=trWtop(:,:,2:length(mygrid.RC));	%3d
    trWbot(:,:,length(mygrid.RC))=0;
    %
    budgO.fluxes.trWtop=trWtop.*RAC3;		%3d
    budgO.fluxes.trWbot=trWbot.*RAC3;		%3d,kg/s
  else
    budgO.fluxes.trWtop=-RAC.*budgO.zconv;	%what does this mean for layer 1 in 3d case?
    budgO.fluxes.trWbot=RAC*0;			%kg/s
  end;

%atn:
  budgO.fluxes.zconv=budgO.fluxes.trWtop-budgO.fluxes.trWbot;	%kg/s , top minus bottom
%done atn

  budgI.fluxes.trWtop=-RAC.*(budgI.zconv+budgO.zconv);	%2d, SIatmFW*RAC [kg/m2/s*m2=kg/s]
  budgI.fluxes.trWbot=-RAC.*budgO.zconv;		%2d, kg/s
  
  budgO.zconv=mk3D_mod(RAC,budgO.zconv).*budgO.zconv;	%2d, kg/s , shouldn't need mk3D
  budgI.zconv=RAC.*budgI.zconv;%kg/s			%2d
  budgOI.zconv=budgO.zconv+budgI.zconv;	%this is what GAel has for top to bottom
  
  tmp=(SIatmFW-oceFWflx-SIsnPrcp).*RAC;	%F(A->I)+F(A->O)-F(I->O)-F(A->O)-F(A->I) = -F(I->O)
  budgO.fluxes.ice2ocean=-tmp;		%2d, kg/s, vol flux ice to ocean, note minus sign
  tmp=(SIatmFW-SIsnPrcp).*RAC;		%F(A->I) + F(A->O) - F(A->I) = F(A->O)
  budgO.fluxes.atm2ocean=tmp;		%2d, kg/s , should be equal to sum(budgO.fluxes.zconv,2)-budgO.fluxes.ice2ocean;
  %confirm: sum(sum(ckO.zconv,2)-(ckO.atm2ocean+ckO.ice2ocean)) = -1.6e-9, yes!

%===================== BLOCK 3: Horizontal convergence =========================================
%ocean:
  varName={'UVELMASS','VVELMASS'};
  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

  if(test3d);
    %3D UVELMASS,VVELMASS are multiplied by DRF /*{{{*/
    %(2D diagnostics are expectedly vertically integrated by MITgcm) %/*}}}*/
    tmp1=drf3d;%mk3D_mod(mygrid.DRF,UVELMASS);
    UVELMASS=tmp1.*UVELMASS;			%3d, m2/s
    VVELMASS=tmp1.*VVELMASS;			%3d, m2/s
    dxg1=dxg3d; dyg1=dyg3d;
  else;
    dxg1=dxg;dyg1=dyg;
  end;

%note: check, before, i had d[y,x]g instead of d[y,x]g1 for the 2 lines below.
  tmpUo=myparms.rhoconst*dyg1.*UVELMASS;		%2d or 3d, kg/m3*m*m2/s = [kg/s]
  tmpVo=myparms.rhoconst*dxg1.*VVELMASS;		%2d or 3d

  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));

%add bolus, only real4.  If i'm correct, [u,v]velmass includes already bolus[u,v]
  if(test3d);
    clear gm_psi* bolus* tmpUo tmpVo tmp 

    if(read_netcdf);
      varName={'GM_PsiX','GM_PsiY'};n3=nz;
      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([lower(strtrim(varName{ifld})) '=squeeze(fldOut(:,:,:,1));']); 
      end;
      clear fldOut
    else;
      gm_psix=read_slice([dirIn '../TRSP/trsp_3d_set2.' sprintf('%10.10i',t2) '.data'],nx,ny,[1:nz]+1*nz);
      gm_psiy=read_slice([dirIn '../TRSP/trsp_3d_set2.' sprintf('%10.10i',t2) '.data'],nx,ny,[1:nz]+2*nz);
    end;
    [bolusU,bolusV,bolusW]=calc_bolus_mod(gm_psix,gm_psiy,nfx,nfy);	%outputs in compact format
    tmpUo=myparms.rhoconst*dyg3d.*drf3d.*bolusU;			%kg/m3*m*m*m/s = [kg/s]
    tmpVo=myparms.rhoconst*dxg3d.*drf3d.*bolusV;			%
    tmp=calc_UV_conv_mod(tmpUo,tmpVo,nfx,nfy);
    budgO.fluxes.hbolus=cat(2,tmp{1},tmp{3},reshape(tmp{4},nfy(4),nfx(4),nz),reshape(tmp{5},nfy(5),nfx(5),nz));
    budgO.fluxes.zbolus=myparms.rhoconst.*bolusW.*RAC3;	%kg/m3*m/s*m2
  end;

  if(test3d)
    tmpUo=myparms.rhoconst*dyg3d.*UVELMASS;		%kg/m3*m*m2/s = [kg/s], note UVELMASS is now UVELMASS*drF
    tmpVo=myparms.rhoconst*dxg3d.*VVELMASS;		%
  end;

%ice
%note for ASTE R1, all the diffusive terms for sea ice are zeros, so we skip reading
  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

  varName={'DFxEHEFF','DFyEHEFF','DFxESNOW','DFyESNOW'};n3=1;
  for ifld=1:size(varName,2);
    fldOut=zeros(nx,ny,n3);eval([varName{ifld} '=fldOut;']);clear fldOut;
  end;
%/*{{{*/
%ADVxHEFF=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,1,'real*8');
%ADVyHEFF=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,2,'real*8');
%DFxEHEFF=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,3,'real*8');
%DFyEHEFF=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,4,'real*8');
%ADVxSNOW=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,5,'real*8');
%ADVySNOW=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,6,'real*8');
%DFxESNOW=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,7,'real*8');
%DFyESNOW=read_slice([dirIn 'budg2d_hflux_set1.' t2str '.data'],nx,ny,8,'real*8');/*}}}*/

  tmpUi=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+...	%m.m2/s * kg/m3 = [kg/s]
         myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);		%
  tmpVi=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+...
         myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);

  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, [kg/s]

  if(test3d)
%this is the right place to compute transports across sections:
    tmpUi=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+...		%m.m2/s * kg/m3 = [kg/s]
           myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);		%
    tmpVi=(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 kg/s

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

  %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_vol_ocn glo_vol_tot glo_vol_ice glo_bp'];	%/*{{{*/ [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_vol_tot_lev';
    end;							%/*}}}*/ [2]
%--------------------

%now do global sum 
%after this sum, unit is: [kg/s/m2]				/*{{{*/ [2]
    tmp1=calc_mskmean_T_mod(budgOI,hfC1p,RACgp,'extensive');glo_vol_tot=[tmp1.tend;tmp1.hconv;tmp1.zconv];
    tmp1=calc_mskmean_T_mod(budgO,hfC1p,RACgp,'extensive'); glo_vol_ocn=[tmp1.tend;tmp1.hconv;tmp1.zconv];
    tmp1=calc_mskmean_T_mod(budgI,hfC1p,RACgp,'extensive'); glo_vol_ice=[tmp1.tend;tmp1.hconv;tmp1.zconv];
    glo_bp=nansum(bp(:).*hfC1p(:).*RACgp(:))/nansum(hfC1p(:).*RACgp(:));

%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_vol_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.01*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.01*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.01*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.01*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
%/*}}}*/
    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;				%/*}}}*/
%============================================
%MISSION ACCOMPLISHED!!! onto heat!
%===================
end;%ts loop
fprintf('\n');
