matlab_class/gcmfaces_diags/diags_set_D.m


%select kBudget:
if ~isempty(setDiagsParams);
  kBudget=setDiagsParams{1};
else;
  kBudget=1;
end;

doMoreBudgetOutput=0;
%doMoreBudgetOutput=1;

%override default file name:
%---------------------------
tmp1=setDiags;
if kBudget>1;
    tmp1=sprintf('D%02i',kBudget);
end;
fileMat=['diags_set_' tmp1];

if userStep==1;%diags to be computed
    listDiags=['glo_vol_ocn glo_vol_tot glo_vol_ice glo_bp'];
    listDiags=[listDiags ' north_vol_ocn north_vol_tot north_vol_ice north_bp'];
    listDiags=[listDiags ' south_vol_ocn south_vol_tot south_vol_ice south_bp'];
    listDiags=[listDiags ' glo_heat_ocn glo_heat_tot glo_heat_ice'];
    listDiags=[listDiags ' north_heat_ocn north_heat_tot north_heat_ice'];
    listDiags=[listDiags ' south_heat_ocn south_heat_tot south_heat_ice'];
    listDiags=[listDiags ' glo_salt_ocn glo_salt_tot glo_salt_ice'];
    listDiags=[listDiags ' north_salt_ocn north_salt_tot north_salt_ice'];
    listDiags=[listDiags ' south_salt_ocn south_salt_tot south_salt_ice'];

elseif userStep==2;%input files and variables
    dirSnap=fullfile(dirModel,'diags',filesep,'BUDG',filesep);
    if ~isdir(dirSnap); dirSnap=fullfile(dirModel,'diags',filesep); end;
    tmp1=fullfile(dirSnap,'budg2d_snap_set2*meta');
    test3d=isempty(dir(tmp1));
    %
    listFlds={    'ETAN','SIheff','SIhsnow','THETA   ','SALT    ','PHIBOT','geothFlux'};
    listFlds={listFlds{:},'SIatmFW ','oceFWflx','SItflux','TFLUX','SFLUX','oceSPflx','SRELAX'};
    listFlds={listFlds{:},'oceQnet ','SIatmQnt','SIaaflux','SIsnPrcp','SIacSubl'};
    listFlds={listFlds{:},'TRELAX','WTHMASS','WSLTMASS','oceSflux','oceQsw','oceSPtnd'};
    if kBudget>1|test3d;
        listFlds={listFlds{:},'ADVr_TH','DFrE_TH','DFrI_TH','ADVr_SLT','DFrE_SLT','DFrI_SLT','WVELMASS'};
    end;
    listFlds={listFlds{:},'SDIAG1','SDIAG2','SDIAG3'};
    listFlds={listFlds{:},'UVELMASS','VVELMASS','AB_gT','AB_gS'};
    listFlds={listFlds{:},'ADVx_TH ','ADVy_TH ','DFxE_TH ','DFyE_TH '};
    listFlds={listFlds{:},'ADVx_SLT','ADVy_SLT','DFxE_SLT','DFyE_SLT'};
    listFlds={listFlds{:},'ADVxHEFF','ADVyHEFF','DFxEHEFF','DFyEHEFF'};
    listFlds={listFlds{:},'ADVxSNOW','ADVySNOW','DFxESNOW','DFyESNOW'};
    listFldsNames=deblank(listFlds);
    %
    listFiles={'rate_budg2d_snap_set1','budg2d_hflux_set1','budg2d_zflux_set1','budg2d_zflux_set2'};
    if test3d;
        listFiles={listFiles{:},'rate_budg3d_snap_set1','budg3d_hflux_set1','budg3d_zflux_set1'};
    elseif kBudget==1;
        listFiles={listFiles{:},'rate_budg2d_snap_set2','budg2d_hflux_set2','geothermalFlux'};
    else;
        tmp1=sprintf('rate_budg2d_snap_set3_%02i',kBudget);
        tmp2=sprintf('budg2d_zflux_set3_%02i',kBudget);
        tmp3=sprintf('budg2d_hflux_set3_%02i',kBudget);
        tmp4=sprintf('geothermalFlux_%02i',kBudget);
        listFiles={listFiles{:},tmp1,tmp2,tmp3,tmp4};
    end;
    listSubdirs={[dirMat 'BUDG/' ],[dirMat '../BUDG/' ],dirSnap};

elseif userStep==3;%computational part;

    %preliminary tests
    test1=isempty(dir([dirModel 'diags/BUDG/budg2d_snap_set1*']));
    test2=isempty(dir([dirMat 'BUDG/rate_budg2d_snap_set1*']))&...
          isempty(dir([dirMat '../BUDG/rate_budg2d_snap_set1*']));
    if (strcmp(setDiags,'D')&test1&test2);
        fprintf('\n abort : global and regional budgets, due to missing \n');
        fprintf(['\n   ' dirModel 'diags/BUDG/budg2d_snap_set1* \n']);
        return;
    end;

    if (strcmp(setDiags,'D')&test2);
        fprintf('\n abort : global and regional budgets, due to missing \n');
        fprintf(['\n   ' dirModel 'diags/BUDG/rate_budg2d_snap_set1* \n']);
        return;
    end;
    
    %override default file name:
    %---------------------------
    tmp1=setDiags;
    if kBudget>1;
        tmp1=sprintf('D%02i',kBudget);
    end;
    fileMat=['diags_set_' tmp1 '_' num2str(tt) '.mat'];
        
    %fill in optional fields:
    %------------------------
    if isempty(who('TRELAX')); TRELAX=0; end;
    if isempty(who('SRELAX')); SRELAX=0; end;
    if isempty(who('AB_gT')); AB_gT=0; end;
    if isempty(who('AB_gS')); AB_gS=0; end;
    if isempty(who('oceSPtnd')); oceSPtnd=0; end;
    if isempty(who('oceSPflx')); oceSPflx=0; end;
    if isempty(who('PHIBOT')); PHIBOT=0; end;
    if isempty(who('geothFlux')); geothFlux=0; end;

    %aliases from development phase (applies to 2012 core runs)
    %---------------------------------------------------------
    if ~isempty(who('SDIAG1')); SRELAX=SDIAG1; end;
    if ~isempty(who('SDIAG2')); SIatmFW=SDIAG2; end;
    if ~isempty(who('SDIAG3')); SItflux=SDIAG3; end;

    %=======indexing and sign convention======

    %- MITgcm: fluxes are >0 downward, k=1 start at free surface
    %- here, similarly: >0 downward, k=1 free surface k=2 sea floor        
    if ~test3d;
      budgO.specs.top='free surface';
      if kBudget>1; budgO.specs.top=['interface no. ' num2str(kBudget)]; end;
      budgO.specs.bottom='sea floor';
    else;
      budgO.specs.top='interface k';
      budgO.specs.bottom='interface k+1';
    end;
    budgI.specs.top='ocn-ice to atm interface';
    budgI.specs.bottom='free surface';
    
    %here we output tendencies and fluxes in kg/s
    budgMo=budgO; budgMi=budgI;
    budgMo.specs.units='kg/s';%ocean only
    budgMi.specs.units='kg/s';%ice only
    %here we output tendencies and fluxes in Watts
    budgHo=budgO; budgHi=budgI;
    budgHo.specs.units='W';%ocean only
    budgHi.specs.units='W';%ice only
    %here we output tendencies and fluxes in g/s
    budgSo=budgO; budgSi=budgI;
    budgSo.specs.units='g/s';%ocean only
    budgSi.specs.units='g/s';%ice only        

    %=======MASS=========
    
    [contOCN,hdivOCN,zdivOCN,budgMo,contICE,hdivICE,zdivICE,budgMi]=...
       calc_budget_mass(kBudget);

    contTOT=contOCN+contICE;
    hdivTOT=hdivOCN+hdivICE;
    zdivTOT=zdivOCN+zdivICE;
    
    %bottom pressure for comparison:
    bp=myparms.rhoconst/9.81*PHIBOT;
    
    %compute global integrals:
    %-------------------------
    msk=mygrid.mskC(:,:,kBudget);
    glo_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    glo_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    glo_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
    glo_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
    
    %compute northern hemisphere integrals:
    msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
    north_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    north_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    north_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
    north_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
    
    %and southern hemisphere integrals:
    msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
    south_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    south_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    south_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
    south_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
    
    %=======HEAT=======

    [contOCN,hdivOCN,zdivOCN,budgMo,contICE,hdivICE,zdivICE,budgMi]=...
       calc_budget_heat(kBudget);

    contTOT=contOCN+contICE;
    hdivTOT=hdivOCN+hdivICE;
    zdivTOT=zdivOCN+zdivICE;
    
    %compute global integrals:
    %-------------------------
    msk=mygrid.mskC(:,:,kBudget);
    glo_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    glo_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    glo_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
    
    %compute northern hemisphere integrals:
    msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
    north_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    north_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    north_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
    
    %and southern hemisphere integrals:
    msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
    south_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    south_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    south_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);

    %=======SALT=======

    [contOCN,hdivOCN,zdivOCN,budgMo,contICE,hdivICE,zdivICE,budgMi]=...
       calc_budget_salt(kBudget);

    contTOT=contOCN+contICE;
    hdivTOT=hdivOCN+hdivICE;
    zdivTOT=zdivOCN+zdivICE;
        
    %compute global integrals:
    %-------------------------
    msk=mygrid.mskC(:,:,kBudget);
    glo_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    glo_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    glo_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
    
    %compute northern hemisphere integrals:
    msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
    north_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    north_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    north_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
    
    %and southern hemisphere integrals:
    msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
    south_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
    south_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
    south_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);

    if doMoreBudgetOutput;
        %list of budgets to output
        listbudg={'budgMo','budgHo','budgSo'};
        if kBudget==1; listbudg={listbudg{:},'budgMi','budgHi','budgSi'}; end;
        %the actual output
        for iibudg=1:length(listbudg);
            %set directory name
            dirbudg=dirMat;
            if ~isempty(strfind(dirMat,['diags_set_' setDiags '/']))
                dirbudg=fullfile(dirMat,'..',filesep);
            end;
            sufbudg=''; 
            if kBudget>1; sufbudg=num2str(kBudget); end;
            dirbudg=fullfile(dirbudg,['diags_set_' listbudg{iibudg} sufbudg],filesep);
            %
            if ~isdir(dirbudg); mkdir(dirbudg); end;
            %set file name
            filebudg=[listbudg{iibudg} '_' num2str(tt) '.mat'];
            %output to file
            eval(['tmpbudg=' listbudg{iibudg} ';']);
            save([dirbudg filebudg],'-struct','tmpbudg');
        end;
    end;

%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
%===================== PLOTTING SEQUENCE BEGINS    =========================%

elseif userStep==-1;%plotting

    if isempty(setDiagsParams);
      choicePlot={'all'};
    elseif isnumeric(setDiagsParams{1})&length(setDiagsParams)==1;
      choicePlot={'all'};
    elseif isnumeric(setDiagsParams{1});
      choicePlot={setDiagsParams{2:end}};
    else;
      choicePlot=setDiagsParams;
    end;

    tt=[1:length(alldiag.listTimes)];
    TT=alldiag.listTimes(tt);
    nt=length(TT);

    if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));

        %1.1) ocean+seaice mass budgets
        %------------------------------
        figureL;
        %global volume budget:
        subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_tot,'kg/m^2','Global Mean Mass (incl. ice)');
        %add bp:
        dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
        plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
        %northern hemisphere budget:
        subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_tot,'kg/m^2','Northern Mean Mass (incl. ice)');
        %add bp:
        dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
        plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
        %southern hemisphere budget:
        subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_tot,'kg/m^2','Southern Mean Mass (incl. ice)');
        %add bp:
        dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
        plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
        %add to tex file
        myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
        myCaption={myCaption{:},'mass budget (ocean+ice) in kg/m$^2$.'};
        if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

        %1.2) ice mass budgets
        %---------------------
        figureL;
        subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ice,'kg/m^2','Global Mean Mass (only ice)');
        dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
        plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
        subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ice,'kg/m^2','Northern Mean Mass (only ice)');
        dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
        plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
        subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ice,'kg/m^2','Southern Mean Mass (only ice)');
        dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
        plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
        %add to tex file
        myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
        myCaption={myCaption{:},'mass budget (ice only) in kg/m$^2$.'};
        if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

    end;

    if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));

    %1.3) ocean mass budgets
    %-----------------------
    figureL;
    %global volume budget:
    subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ocn,'kg/m^2','Global Mean Mass (only ocean)');
    dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
    plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
    subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ocn,'kg/m^2','Northern Mean Mass (only ocean)');
    dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
    plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
    subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ocn,'kg/m^2','Southern Mean Mass (only ocean)');
    dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
    plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
    %add to tex file
    myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
    myCaption={myCaption{:},'mass budget (ocean only) in kg/m$^2$.'};
    if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

    end;

    if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));

        %2.1) ocean+seaice heat budgets
        %------------------------------
        figureL;
        subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_tot,'J/m^2','Global Mean Ocean Heat (incl. ice)');
        subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_tot,'J/m^2','Northern Mean Ocean Heat (incl. ice)');
        subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_tot,'J/m^2','Southern Mean Ocean Heat (incl. ice)');
        %add to tex file
        myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
        myCaption={myCaption{:},'heat budget (ocean+ice) in J/m$^2$.'};
        if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

        %2.2) ice heat budgets
        %---------------------
        figureL;
        subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_ice,'J/m^2','Global Mean Ocean Heat (only ice)');
        subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ice,'J/m^2','Northern Mean Ocean Heat (only ice)');
        subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ice,'J/m^2','Southern Mean Ocean Heat (only ice)');
        %add to tex file
        myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
        myCaption={myCaption{:},'heat budget (ice only) in J/m$^2$.'};
        if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

    end;

    if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));

    %2.3) ocean heat budgets
    %-----------------------
    figureL;
    subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_ocn,'J/m^2','Global Mean Ocean Heat (only ocean)');
    subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ocn,'J/m^2','Northern Mean Ocean Heat (only ocean)');
    subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ocn,'J/m^2','Southern Mean Ocean Heat (only ocean)');
    %add to tex file
    myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
    myCaption={myCaption{:},'heat budget (ocean only) in J/m$^2$.'};
    if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

    end;

    if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));

        %3.1) ocean+seaice salt budgets
        %------------------------------
        figureL;
        subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_tot,'g/m^2','Global Mean Ocean Salt (incl. ice)');
        subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_tot,'g/m^2','Northern Mean Ocean Salt (incl. ice)');
        subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_tot,'g/m^2','Southern Mean Ocean Salt (incl. ice)');
        %add to tex file
        myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
        myCaption={myCaption{:},'salt budget (ocean+ice) in g/m$^2$.'};
        if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

        %2.2) ice salt budgets
        %---------------------
        figureL;
        subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_ice,'g/m^2','Global Mean Ocean Salt (only ice)');
        subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ice,'g/m^2','Northern Mean Ocean Salt (only ice)');
        subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ice,'g/m^2','Southern Mean Ocean Salt (only ice)');
        %add to tex file
        myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
        myCaption={myCaption{:},'salt budget (ice only) in g/m$^2$.'};
        if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

    end;


    if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));

    %3.3) ocean salt budgets
    %-----------------------
    figureL;
    subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_ocn,'g/m^2','Global Mean Ocean Salt (only ocean)');
    subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ocn,'g/m^2','Northern Mean Ocean Salt (only ocean)');
    subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ocn,'g/m^2','Southern Mean Ocean Salt (only ocean)');
    %add to tex file
    myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
    myCaption={myCaption{:},'salt budget (ocean only) in g/m$^2$.'};
    if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;

    end;

end;
1	gforget	1.1
2			%select kBudget:
3			if ~isempty(setDiagsParams);
4			kBudget=setDiagsParams{1};
5			else;
6			kBudget=1;
7			end;
8
9	gforget	1.9	doMoreBudgetOutput=0;
10	gforget	1.16	%doMoreBudgetOutput=1;
11	gforget	1.9
12	gforget	1.2	%override default file name:
13			%---------------------------
14			tmp1=setDiags;
15			if kBudget>1;
16			tmp1=sprintf('D%02i',kBudget);
17	gforget	1.1	end;
18	gforget	1.2	fileMat=['diags_set_' tmp1];
19	gforget	1.1
20			if userStep==1;%diags to be computed
21			listDiags=['glo_vol_ocn glo_vol_tot glo_vol_ice glo_bp'];
22			listDiags=[listDiags ' north_vol_ocn north_vol_tot north_vol_ice north_bp'];
23			listDiags=[listDiags ' south_vol_ocn south_vol_tot south_vol_ice south_bp'];
24			listDiags=[listDiags ' glo_heat_ocn glo_heat_tot glo_heat_ice'];
25			listDiags=[listDiags ' north_heat_ocn north_heat_tot north_heat_ice'];
26			listDiags=[listDiags ' south_heat_ocn south_heat_tot south_heat_ice'];
27			listDiags=[listDiags ' glo_salt_ocn glo_salt_tot glo_salt_ice'];
28			listDiags=[listDiags ' north_salt_ocn north_salt_tot north_salt_ice'];
29			listDiags=[listDiags ' south_salt_ocn south_salt_tot south_salt_ice'];
30	gforget	1.2
31	gforget	1.1	elseif userStep==2;%input files and variables
32	gforget	1.15	dirSnap=fullfile(dirModel,'diags',filesep,'BUDG',filesep);
33			if ~isdir(dirSnap); dirSnap=fullfile(dirModel,'diags',filesep); end;
34			tmp1=fullfile(dirSnap,'budg2d_snap_set2*meta');
35	gforget	1.11	test3d=isempty(dir(tmp1));
36			%
37	gforget	1.16	listFlds={ 'ETAN','SIheff','SIhsnow','THETA ','SALT ','PHIBOT','geothFlux'};
38	gforget	1.1	listFlds={listFlds{:},'SIatmFW ','oceFWflx','SItflux','TFLUX','SFLUX','oceSPflx','SRELAX'};
39			listFlds={listFlds{:},'oceQnet ','SIatmQnt','SIaaflux','SIsnPrcp','SIacSubl'};
40			listFlds={listFlds{:},'TRELAX','WTHMASS','WSLTMASS','oceSflux','oceQsw','oceSPtnd'};
41	gforget	1.11	if kBudget>1\|test3d;
42	gforget	1.1	listFlds={listFlds{:},'ADVr_TH','DFrE_TH','DFrI_TH','ADVr_SLT','DFrE_SLT','DFrI_SLT','WVELMASS'};
43			end;
44	gforget	1.4	listFlds={listFlds{:},'SDIAG1','SDIAG2','SDIAG3'};
45	gforget	1.1	listFlds={listFlds{:},'UVELMASS','VVELMASS','AB_gT','AB_gS'};
46			listFlds={listFlds{:},'ADVx_TH ','ADVy_TH ','DFxE_TH ','DFyE_TH '};
47			listFlds={listFlds{:},'ADVx_SLT','ADVy_SLT','DFxE_SLT','DFyE_SLT'};
48			listFlds={listFlds{:},'ADVxHEFF','ADVyHEFF','DFxEHEFF','DFyEHEFF'};
49			listFlds={listFlds{:},'ADVxSNOW','ADVySNOW','DFxESNOW','DFyESNOW'};
50			listFldsNames=deblank(listFlds);
51			%
52			listFiles={'rate_budg2d_snap_set1','budg2d_hflux_set1','budg2d_zflux_set1','budg2d_zflux_set2'};
53	gforget	1.11	if test3d;
54			listFiles={listFiles{:},'rate_budg3d_snap_set1','budg3d_hflux_set1','budg3d_zflux_set1'};
55			elseif kBudget==1;
56	gforget	1.16	listFiles={listFiles{:},'rate_budg2d_snap_set2','budg2d_hflux_set2','geothermalFlux'};
57	gforget	1.1	else;
58			tmp1=sprintf('rate_budg2d_snap_set3_%02i',kBudget);
59			tmp2=sprintf('budg2d_zflux_set3_%02i',kBudget);
60			tmp3=sprintf('budg2d_hflux_set3_%02i',kBudget);
61	gforget	1.16	tmp4=sprintf('geothermalFlux_%02i',kBudget);
62			listFiles={listFiles{:},tmp1,tmp2,tmp3,tmp4};
63	gforget	1.1	end;
64	gforget	1.15	listSubdirs={[dirMat 'BUDG/' ],[dirMat '../BUDG/' ],dirSnap};
65	gforget	1.2
66	gforget	1.1	elseif userStep==3;%computational part;
67	gforget	1.2
68			%preliminary tests
69			test1=isempty(dir([dirModel 'diags/BUDG/budg2d_snap_set1*']));
70	gforget	1.5	test2=isempty(dir([dirMat 'BUDG/rate_budg2d_snap_set1*']))&...
71			isempty(dir([dirMat '../BUDG/rate_budg2d_snap_set1*']));
72	gforget	1.2	if (strcmp(setDiags,'D')&test1&test2);
73			fprintf('\n abort : global and regional budgets, due to missing \n');
74			fprintf(['\n ' dirModel 'diags/BUDG/budg2d_snap_set1* \n']);
75			return;
76			end;
77
78			if (strcmp(setDiags,'D')&test2);
79			fprintf('\n abort : global and regional budgets, due to missing \n');
80			fprintf(['\n ' dirModel 'diags/BUDG/rate_budg2d_snap_set1* \n']);
81			return;
82			end;
83	gforget	1.1
84			%override default file name:
85			%---------------------------
86			tmp1=setDiags;
87			if kBudget>1;
88			tmp1=sprintf('D%02i',kBudget);
89			end;
90			fileMat=['diags_set_' tmp1 '_' num2str(tt) '.mat'];
91
92			%fill in optional fields:
93			%------------------------
94	gforget	1.11	if isempty(who('TRELAX')); TRELAX=0; end;
95			if isempty(who('SRELAX')); SRELAX=0; end;
96			if isempty(who('AB_gT')); AB_gT=0; end;
97			if isempty(who('AB_gS')); AB_gS=0; end;
98			if isempty(who('oceSPtnd')); oceSPtnd=0; end;
99			if isempty(who('oceSPflx')); oceSPflx=0; end;
100			if isempty(who('PHIBOT')); PHIBOT=0; end;
101	gforget	1.16	if isempty(who('geothFlux')); geothFlux=0; end;
102	gforget	1.4
103			%aliases from development phase (applies to 2012 core runs)
104			%---------------------------------------------------------
105			if ~isempty(who('SDIAG1')); SRELAX=SDIAG1; end;
106			if ~isempty(who('SDIAG2')); SIatmFW=SDIAG2; end;
107			if ~isempty(who('SDIAG3')); SItflux=SDIAG3; end;
108
109	gforget	1.17	%=======indexing and sign convention======
110
111			%- MITgcm: fluxes are >0 downward, k=1 start at free surface
112			%- here, similarly: >0 downward, k=1 free surface k=2 sea floor
113			if ~test3d;
114			budgO.specs.top='free surface';
115			if kBudget>1; budgO.specs.top=['interface no. ' num2str(kBudget)]; end;
116			budgO.specs.bottom='sea floor';
117			else;
118			budgO.specs.top='interface k';
119			budgO.specs.bottom='interface k+1';
120			end;
121			budgI.specs.top='ocn-ice to atm interface';
122			budgI.specs.bottom='free surface';
123	gforget	1.1
124	gforget	1.17	%here we output tendencies and fluxes in kg/s
125			budgMo=budgO; budgMi=budgI;
126			budgMo.specs.units='kg/s';%ocean only
127			budgMi.specs.units='kg/s';%ice only
128			%here we output tendencies and fluxes in Watts
129			budgHo=budgO; budgHi=budgI;
130			budgHo.specs.units='W';%ocean only
131			budgHi.specs.units='W';%ice only
132			%here we output tendencies and fluxes in g/s
133			budgSo=budgO; budgSi=budgI;
134			budgSo.specs.units='g/s';%ocean only
135			budgSi.specs.units='g/s';%ice only
136
137	gforget	1.1	%=======MASS=========
138
139	gforget	1.17	[contOCN,hdivOCN,zdivOCN,budgMo,contICE,hdivICE,zdivICE,budgMi]=...
140			calc_budget_mass(kBudget);
141
142	gforget	1.1	contTOT=contOCN+contICE;
143	gforget	1.17	hdivTOT=hdivOCN+hdivICE;
144	gforget	1.1	zdivTOT=zdivOCN+zdivICE;
145
146			%bottom pressure for comparison:
147			bp=myparms.rhoconst/9.81*PHIBOT;
148
149			%compute global integrals:
150			%-------------------------
151			msk=mygrid.mskC(:,:,kBudget);
152			glo_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
153			glo_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
154			glo_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
155			glo_bp=nansum(bp.msk.mygrid.RAC)/nansum(msk.*mygrid.RAC);
156
157			%compute northern hemisphere integrals:
158			msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
159			north_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
160			north_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
161			north_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
162			north_bp=nansum(bp.msk.mygrid.RAC)/nansum(msk.*mygrid.RAC);
163
164			%and southern hemisphere integrals:
165			msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
166			south_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
167			south_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
168			south_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
169			south_bp=nansum(bp.msk.mygrid.RAC)/nansum(msk.*mygrid.RAC);
170
171			%=======HEAT=======
172	gforget	1.11
173	gforget	1.17	[contOCN,hdivOCN,zdivOCN,budgMo,contICE,hdivICE,zdivICE,budgMi]=...
174			calc_budget_heat(kBudget);
175
176	gforget	1.11	contTOT=contOCN+contICE;
177	gforget	1.17	hdivTOT=hdivOCN+hdivICE;
178	gforget	1.1	zdivTOT=zdivOCN+zdivICE;
179
180			%compute global integrals:
181			%-------------------------
182			msk=mygrid.mskC(:,:,kBudget);
183			glo_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
184			glo_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
185			glo_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
186
187			%compute northern hemisphere integrals:
188			msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
189			north_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
190			north_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
191			north_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
192
193			%and southern hemisphere integrals:
194			msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
195			south_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
196			south_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
197			south_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
198	gforget	1.11
199	gforget	1.1	%=======SALT=======
200	gforget	1.17
201			[contOCN,hdivOCN,zdivOCN,budgMo,contICE,hdivICE,zdivICE,budgMi]=...
202			calc_budget_salt(kBudget);
203
204	gforget	1.11	contTOT=contOCN+contICE;
205	gforget	1.17	hdivTOT=hdivOCN+hdivICE;
206	gforget	1.1	zdivTOT=zdivOCN+zdivICE;
207	gforget	1.17
208	gforget	1.1	%compute global integrals:
209			%-------------------------
210			msk=mygrid.mskC(:,:,kBudget);
211			glo_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
212			glo_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
213			glo_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
214
215			%compute northern hemisphere integrals:
216			msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
217			north_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
218			north_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
219			north_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
220
221			%and southern hemisphere integrals:
222			msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
223			south_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
224			south_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
225			south_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
226	gforget	1.2
227	gforget	1.9	if doMoreBudgetOutput;
228			%list of budgets to output
229			listbudg={'budgMo','budgHo','budgSo'};
230			if kBudget==1; listbudg={listbudg{:},'budgMi','budgHi','budgSi'}; end;
231			%the actual output
232			for iibudg=1:length(listbudg);
233			%set directory name
234			dirbudg=dirMat;
235			if ~isempty(strfind(dirMat,['diags_set_' setDiags '/']))
236			dirbudg=fullfile(dirMat,'..',filesep);
237			end;
238			sufbudg='';
239			if kBudget>1; sufbudg=num2str(kBudget); end;
240			dirbudg=fullfile(dirbudg,['diags_set_' listbudg{iibudg} sufbudg],filesep);
241			%
242			if ~isdir(dirbudg); mkdir(dirbudg); end;
243			%set file name
244			filebudg=[listbudg{iibudg} '_' num2str(tt) '.mat'];
245			%output to file
246			eval(['tmpbudg=' listbudg{iibudg} ';']);
247			save([dirbudg filebudg],'-struct','tmpbudg');
248			end;
249			end;
250	gforget	1.13
251	gforget	1.2	%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
252			%===================== PLOTTING SEQUENCE BEGINS =========================%
253
254			elseif userStep==-1;%plotting
255
256			if isempty(setDiagsParams);
257			choicePlot={'all'};
258			elseif isnumeric(setDiagsParams{1})&length(setDiagsParams)==1;
259			choicePlot={'all'};
260			elseif isnumeric(setDiagsParams{1});
261			choicePlot={setDiagsParams{2:end}};
262			else;
263			choicePlot=setDiagsParams;
264			end;
265
266	gforget	1.3	tt=[1:length(alldiag.listTimes)];
267			TT=alldiag.listTimes(tt);
268			nt=length(TT);
269
270	gforget	1.2	if (kBudget==1)&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'mass')));
271
272			%1.1) ocean+seaice mass budgets
273			%------------------------------
274			figureL;
275			%global volume budget:
276	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_tot,'kg/m^2','Global Mean Mass (incl. ice)');
277	gforget	1.2	%add bp:
278			dt=median(diff(TT))86400; bp=dtcumsum(alldiag.glo_bp);
279			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
280			%northern hemisphere budget:
281	heimbach	1.6	subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_tot,'kg/m^2','Northern Mean Mass (incl. ice)');
282	gforget	1.2	%add bp:
283			dt=median(diff(TT))86400; bp=dtcumsum(alldiag.north_bp);
284			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
285			%southern hemisphere budget:
286	heimbach	1.6	subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_tot,'kg/m^2','Southern Mean Mass (incl. ice)');
287	gforget	1.2	%add bp:
288			dt=median(diff(TT))86400; bp=dtcumsum(alldiag.south_bp);
289			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
290			%add to tex file
291			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
292	heimbach	1.6	myCaption={myCaption{:},'mass budget (ocean+ice) in kg/m$^2$.'};
293	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
294
295			%1.2) ice mass budgets
296			%---------------------
297			figureL;
298	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ice,'kg/m^2','Global Mean Mass (only ice)');
299	gforget	1.2	dt=median(diff(TT))86400; bp=dtcumsum(alldiag.glo_bp);
300			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
301	heimbach	1.6	subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ice,'kg/m^2','Northern Mean Mass (only ice)');
302	gforget	1.2	dt=median(diff(TT))86400; bp=dtcumsum(alldiag.north_bp);
303			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
304	heimbach	1.6	subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ice,'kg/m^2','Southern Mean Mass (only ice)');
305	gforget	1.2	dt=median(diff(TT))86400; bp=dtcumsum(alldiag.south_bp);
306			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
307			%add to tex file
308			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
309	heimbach	1.6	myCaption={myCaption{:},'mass budget (ice only) in kg/m$^2$.'};
310	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
311
312			end;
313
314			if (sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'mass')));
315
316			%1.3) ocean mass budgets
317			%-----------------------
318			figureL;
319			%global volume budget:
320	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ocn,'kg/m^2','Global Mean Mass (only ocean)');
321	gforget	1.2	dt=median(diff(TT))86400; bp=dtcumsum(alldiag.glo_bp);
322			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
323	heimbach	1.6	subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ocn,'kg/m^2','Northern Mean Mass (only ocean)');
324	gforget	1.2	dt=median(diff(TT))86400; bp=dtcumsum(alldiag.north_bp);
325			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
326	heimbach	1.6	subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ocn,'kg/m^2','Southern Mean Mass (only ocean)');
327	gforget	1.2	dt=median(diff(TT))86400; bp=dtcumsum(alldiag.south_bp);
328			plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
329			%add to tex file
330			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
331	heimbach	1.6	myCaption={myCaption{:},'mass budget (ocean only) in kg/m$^2$.'};
332	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
333
334			end;
335
336			if (kBudget==1)&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'heat')));
337
338			%2.1) ocean+seaice heat budgets
339			%------------------------------
340			figureL;
341	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_tot,'J/m^2','Global Mean Ocean Heat (incl. ice)');
342			subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_tot,'J/m^2','Northern Mean Ocean Heat (incl. ice)');
343			subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_tot,'J/m^2','Southern Mean Ocean Heat (incl. ice)');
344	gforget	1.2	%add to tex file
345			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
346	heimbach	1.6	myCaption={myCaption{:},'heat budget (ocean+ice) in J/m$^2$.'};
347	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
348
349			%2.2) ice heat budgets
350			%---------------------
351			figureL;
352	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_ice,'J/m^2','Global Mean Ocean Heat (only ice)');
353			subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ice,'J/m^2','Northern Mean Ocean Heat (only ice)');
354			subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ice,'J/m^2','Southern Mean Ocean Heat (only ice)');
355	gforget	1.2	%add to tex file
356			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
357	heimbach	1.6	myCaption={myCaption{:},'heat budget (ice only) in J/m$^2$.'};
358	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
359
360			end;
361
362			if (sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'heat')));
363
364			%2.3) ocean heat budgets
365			%-----------------------
366			figureL;
367	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_ocn,'J/m^2','Global Mean Ocean Heat (only ocean)');
368			subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ocn,'J/m^2','Northern Mean Ocean Heat (only ocean)');
369			subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ocn,'J/m^2','Southern Mean Ocean Heat (only ocean)');
370	gforget	1.2	%add to tex file
371			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
372	heimbach	1.6	myCaption={myCaption{:},'heat budget (ocean only) in J/m$^2$.'};
373	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
374
375			end;
376
377			if (kBudget==1)&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'salt')));
378
379			%3.1) ocean+seaice salt budgets
380			%------------------------------
381			figureL;
382	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_tot,'g/m^2','Global Mean Ocean Salt (incl. ice)');
383			subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_tot,'g/m^2','Northern Mean Ocean Salt (incl. ice)');
384			subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_tot,'g/m^2','Southern Mean Ocean Salt (incl. ice)');
385	gforget	1.2	%add to tex file
386			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
387	heimbach	1.6	myCaption={myCaption{:},'salt budget (ocean+ice) in g/m$^2$.'};
388	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
389
390			%2.2) ice salt budgets
391			%---------------------
392			figureL;
393	heimbach	1.6	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_ice,'g/m^2','Global Mean Ocean Salt (only ice)');
394			subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ice,'g/m^2','Northern Mean Ocean Salt (only ice)');
395			subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ice,'g/m^2','Southern Mean Ocean Salt (only ice)');
396	gforget	1.2	%add to tex file
397			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
398	heimbach	1.6	myCaption={myCaption{:},'salt budget (ice only) in g/m$^2$.'};
399	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
400
401			end;
402
403
404			if (sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'salt')));
405
406			%3.3) ocean salt budgets
407			%-----------------------
408			figureL;
409	gforget	1.8	subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_ocn,'g/m^2','Global Mean Ocean Salt (only ocean)');
410			subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ocn,'g/m^2','Northern Mean Ocean Salt (only ocean)');
411			subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ocn,'g/m^2','Southern Mean Ocean Salt (only ocean)');
412	gforget	1.2	%add to tex file
413			myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
414	heimbach	1.6	myCaption={myCaption{:},'salt budget (ocean only) in g/m$^2$.'};
415	gforget	1.2	if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
416
417			end;
418
419	gforget	1.1	end;