%'budgMo' ocean mass budget (volume*rhoconst)
%'budgHo' ocean heat budget
%'budgSo' ocean salt budget
%'budgMi' seaice+snow mass budget
%'budgHi' seaice+snow heat budget
%'budgSi' seaice+snow salt budget

budgName='budgHo'
nt=238
nk=10;
doRes=0;

dirIn=[dir0 '/release1/nctiles_budg/' budgName '/'];

%constant parameters:
parms.rhoconst =1029; %sea water density
parms.rcp      =3994*parms.rhoconst; % sea water rho X heat capacity
parms.rhoi     = 910; %sea ice density
parms.rhosn    = 330; %snow density
parms.flami    = 3.34e05; % latent heat of fusion of ice/snow (J/kg)
parms.flamb    = 2.50e06; % latent heat of evaporation (J/kg)

%reference grid cell volumes:
tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
tmp3=mk3D(mygrid.RAC,tmp1);
vol0=tmp3.*(tmp1);

switch budgName;
case 'budgMo'; fac=parms.rhoconst;
case 'budgHo'; fac=parms.rcp;
case 'budgSo'; fac=parms.rhoconst;
case 'budgMi'; FACheff=myparms.rhoi; FACsnow=myparms.rhosn;
case 'budgHi'; FACheff=-myparms.flami*myparms.rhoi; FACsnow=-myparms.flami*myparms.rhosn;
case 'budgSi'; FACheff=myparms.SIsal0*myparms.rhoi; FACsnow=0;
end;

budgIn=[];

budgIn.ini=read_nctiles([dirIn 'initial/snapshot'],'snapshot');
budgIn.fin=read_nctiles([dirIn 'final/snapshot'],'snapshot');
eval(['ncload ' dirIn 'tend/tend.0001.nc t0']);
eval(['ncload ' dirIn 'tend/tend.0001.nc t1']);
budgIn.dt=t1-t0;

increments=repmat(NaN*mygrid.RAC,[1 1 nt]);

switch budgName;
case 'budgMo'; 
  ini=nansum(budgIn.ini-fac*vol0,3)./mygrid.RAC/fac;
  fin=nansum(budgIn.fin-fac*vol0,3)./mygrid.RAC/fac;
case 'budgHo';
  ini=nansum(budgIn.ini(:,:,1:nk),3)./nansum(vol0(:,:,1:nk),3)/fac;
  fin=nansum(budgIn.fin(:,:,1:nk),3)./nansum(vol0(:,:,1:nk),3)/fac;
otherwise; error('not implemented yet');
end;

for tt=1:nt;
  disp(tt)
  %load the various fields:
  budgIn.tend=read_nctiles([dirIn 'tend/tend'],'tend',tt);
  budgIn.trU=read_nctiles([dirIn 'trU/trU'],'trU',tt);
  budgIn.trV=read_nctiles([dirIn 'trV/trV'],'trV',tt);
  if dirIn(end-1)=='o';
    budgIn.trW=read_nctiles([dirIn 'trW/trW'],'trW',tt);
    budgIn.trWtop=budgIn.trW;
    budgIn.trWbot=budgIn.trW(:,:,2:50);
    budgIn.trWbot(:,:,50)=0;
  else;
    budgIn.trWtop=read_nctiles([dirIn 'trWtop/trWtop'],'trWtop',tt);
    budgIn.trWbot=read_nctiles([dirIn 'trWbot/trWbot'],'trWbot',tt);
  end;
  %compute budget residuals:
  if doRes;
    nr=size(budgIn.tend{1},3);
    for kk=1:nr; prec(kk,:)=check_budg(budgIn,kk); end;
    if tt==1;
      store_prec=prec;
    else;
      store_prec(:,:,tt)=prec;
    end;
  end;
  %reconstruct time series from budget
  switch budgName;
  case 'budgMo'; 
    increments(:,:,tt)=budgIn.dt(tt)*nansum(budgIn.tend,3)./mygrid.RAC/fac;
  case 'budgHo';
    increments(:,:,tt)=budgIn.dt(tt)*nansum(budgIn.tend(:,:,1:nk),3)...
                       ./nansum(vol0(:,:,1:nk),3)/fac;
  otherwise; error('not implemented yet');
  end;

end;%for tt=[1 119 238];