matlab_class/gcmfaces_diags/diags_set_C.m


if userStep==1;%diags to be computed
    listDiags='fldTzonmean fldSzonmean fldETANzonmean fldSSHzonmean fldSIzonmean fldMLDzonmean';
    listDiags=[listDiags ' IceAreaNorth IceAreaSouth IceVolNorth IceVolSouth SnowVolNorth SnowVolSouth'];
    listDiags=[listDiags ' Ueq Teq SSHglo Tglo Sglo TgloProf SgloProf'];
elseif userStep==2;%input files and variables
    listFlds={    'THETA','SALT','ETAN','SIarea','SIheff','SIhsnow','sIceLoad','MXLDEPTH','UVELMASS','VVELMASS'};
    listFldsNames=deblank(listFlds);
    listFiles={'state_2d_set1','other_2d_set1','state_3d_set1','trsp_3d_set1'};
    listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/'],[dirModel 'diags/']};
elseif userStep==3;%computational part;
    %mask and re-arrange fields:
    fldT=THETA.*mygrid.mskC; fldS=SALT.*mygrid.mskC;
    fldETAN=ETAN.*mygrid.mskC(:,:,1);
    fldSSH=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
    fldSI=SIarea.*mygrid.mskC(:,:,1);
    fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
    ux=UVELMASS.*mygrid.mskW; vy=VVELMASS.*mygrid.mskS;
    [ue,un]=calc_UEVNfromUXVY(ux,vy);
    
    %compute zonal means:
    [fldTzonmean]=calc_zonmean_T(fldT);
    [fldSzonmean]=calc_zonmean_T(fldS);
    [fldETANzonmean]=calc_zonmean_T(fldETAN);
    [fldSSHzonmean]=calc_zonmean_T(fldSSH);
    [fldSIzonmean]=calc_zonmean_T(fldSI);
    [fldMLDzonmean]=calc_zonmean_T(fldMLD);
   
    %compute hemispheric ice sums:
    fld=SIarea.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorth=nansum(fld);
    fld=SIarea.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouth=nansum(fld);
    fld=SIheff.*mygrid.RAC.*(mygrid.YC>0); IceVolNorth=nansum(fld);
    fld=SIheff.*mygrid.RAC.*(mygrid.YC<0); IceVolSouth=nansum(fld);
    fld=SIhsnow.*mygrid.RAC.*(mygrid.YC>0); SnowVolNorth=nansum(fld);
    fld=SIhsnow.*mygrid.RAC.*(mygrid.YC<0); SnowVolSouth=nansum(fld);
    
    %equatorial sections of T and U:
    [secX,secY,Teq]=gcmfaces_section([],0,fldT,1);
    [secX,secY,Ueq]=gcmfaces_section([],0,ue,1);
    
    %global means and profiles:
    msk=mygrid.RAC.*mygrid.mskC(:,:,1);
    SSHglo=nansum(fldSSH.*msk)./nansum(msk);
    %
    msk=mygrid.mskC.*mygrid.hFacC;
    msk=msk.*mk3D(mygrid.RAC,msk).*mk3D(mygrid.DRF,msk);
    Tglo=nansum(fldT.*msk)./nansum(msk);
    Sglo=nansum(fldS.*msk)./nansum(msk);
    nr=length(mygrid.RC);
    TgloProf=zeros(nr,1); SgloProf=zeros(nr,1);
    for kk=1:nr;
        rac=mygrid.RAC.*mygrid.mskC(:,:,kk);
        TgloProf(kk)=nansum(rac.*fldT(:,:,kk))/nansum(rac);
        SgloProf(kk)=nansum(rac.*fldS(:,:,kk))/nansum(rac);
    end;


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

elseif userStep==-1;%plotting

    if isempty(setDiagsParams);
      choicePlot={'all'};
    else;
      choicePlot=setDiagsParams;
    end;

    %ensure backward compatibility
    if ~isfield(alldiag,'fldSSHzonmean'); alldiag.fldSSHzonmean=alldiag.fldETANLEADSzonmean; end;
    if ~isfield(alldiag,'SSHglo'); alldiag.SSHglo=NaN*alldiag.Tglo; end;

    %determine number of years in alldiag.listTimes
    myTimes=alldiag.listTimes; 
    %determine the number of records in one year (lYear)
    tmp1=mean(myTimes(2:end)-myTimes(1:end-1));
    lYear=round(1/tmp1);
    %in case when lYear<2 we use records as years
    if ~(lYear>=2); lYear=1; myTimes=[1:length(myTimes)]; end;
    %determine the number of full years (nYears)
    nYears=floor(length(myTimes)/lYear);

    if nYears>1&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'zmtend')));

    if addToTex; write2tex(fileTex,1,'zonal mean tendencies',2); end;

    figureL; set(gcf,'Renderer','zbuffer');
    %cc=[-2:0.25:2];
    cc=1/100*[[-200:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:200]];
    if doAnomalies; cc=scaleAnom/10*cc; end;
    X=mygrid.LATS*ones(1,length(mygrid.RC)); Y=ones(length(mygrid.LATS),1)*(mygrid.RC');
    %T
    subplot(2,1,1);
    fld=annualmean(alldiag.listTimes,alldiag.fldTzonmean,'last')-annualmean(alldiag.listTimes,alldiag.fldTzonmean,'first');
    depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
    cbar=gcmfaces_cmap_cbar(cc); title('zonal mean T anomaly');
    %S
    subplot(2,1,2);
    fld=annualmean(alldiag.listTimes,alldiag.fldSzonmean,'last')-annualmean(alldiag.listTimes,alldiag.fldSzonmean,'first');
    depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
    cbar=gcmfaces_cmap_cbar(cc/4); title('zonal mean S anomaly');
    %
    myCaption={myYmeanTxt,', last year minus first year -- zonal mean temperature (degC; top) and salinity (psu; bottom)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    end;
   
    if sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'eq'));

    if addToTex; write2tex(fileTex,1,'equatorial sections',2); end;

    %time mean equator sections :
    figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2];
    [secX,secY,LONeq]=gcmfaces_section([],0,mygrid.XC,1);
    X=LONeq*ones(1,length(mygrid.RC)); Y=ones(length(LONeq),1)*mygrid.RC';
    X=circshift(X,[-180 0]); X(X<0)=X(X<0)+360;
    %T
    subplot(2,1,1);
    fld=annualmean(alldiag.listTimes,alldiag.Teq,'all'); fld=circshift(fld,[-180 0]);
    depthStretchPlot('pcolor',{X,Y,fld},[0:50:350],[0 350 350]); shading interp;
    if ~doAnomalies; cc=18+6*cc0; else; cc=scaleAnom/10*cc0*2; end;
    cbar=gcmfaces_cmap_cbar(cc); title('equator T (degC)');
    %U
    subplot(2,1,2);
    fld=annualmean(alldiag.listTimes,alldiag.Ueq,'all'); fld=circshift(fld,[-180 0]);
    depthStretchPlot('pcolor',{X,Y,fld},[0:50:350],[0 350 350]); shading interp;
    if ~doAnomalies; cc=cc0/5*2; else; cc=scaleAnom/10*cc0/5; end;
    cbar=gcmfaces_cmap_cbar(cc); title('equator zonal velocity (m/s)');
    %
    myCaption={myYmeanTxt,' mean -- equator temperature (degC;top) and zonal velocity (m/s;bottom)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    end;

    if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'gmtime')));

    if addToTex; write2tex(fileTex,1,'global mean properties',2); end;

    %global mean T/S:
    figureL;
    subplot(3,1,1); vec=alldiag.SSHglo;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
    aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
    legend('monthly','annual mean'); title('Global Mean Sea level (m, uncorrected free surface)');
    subplot(3,1,2); vec=alldiag.Tglo;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
    aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
    legend('monthly','annual mean'); title('Global Mean Temperature (degC)');
    subplot(3,1,3); vec=alldiag.Sglo;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
    aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
    legend('monthly','annual mean'); title('Global Mean Salinity (psu)');
    myCaption={'global mean T (degC; top) and S (psu; bottom)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
   
    %global mean T/S profiles:
    figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2];
    X=TT*ones(1,length(mygrid.RC)); Y=ones(length(TT),1)*(mygrid.RC'); X=X'; Y=Y';
    %T
    subplot(2,1,1);
    [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.TgloProf),'first');
    depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
    if ~doAnomalies; cc=cc0/5; else; cc=scaleAnom/10*cc0/10; end;
    cbar=gcmfaces_cmap_cbar(cc); title('global mean T minus first year (K)');
    %S
    subplot(2,1,2);
    [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.SgloProf),'first');
    depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
    if ~doAnomalies; cc=cc0/50; else; cc=scaleAnom/10*cc0/50; end;
    cbar=gcmfaces_cmap_cbar(cc); title('global mean S minus first year (psu)');
    %
    myCaption={'global mean temperature (K; top) and salinity (psu; bottom) minus first year'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    end;


    if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'lzmtime')));

    if addToTex; write2tex(fileTex,1,'zonal mean properties',2); end;

    listLatPlot=[-75 -65 -45 -25 0 25 45 65 75];
    %zonal mean temperature profile at selected latitude:
    for iLatPlot=1:length(listLatPlot);
        tmp1=abs(mygrid.LATS-listLatPlot(iLatPlot));
        iLat=find(tmp1==min(tmp1)); iLat=iLat(1);
        figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2];
        if doAnomalies; cc=scaleAnom/10*cc; end;
        %T
        subplot(2,1,1);
        [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.fldTzonmean(iLat,:,:)),'first');
        X=TT*ones(1,length(mygrid.RC)); Y=ones(length(TT),1)*(mygrid.RC'); X=X'; Y=Y';
        title0=['mean T minus first year (K) at lat ~ ' num2str(listLatPlot(iLatPlot))];
        depthStretchPlot('pcolor',{X,Y,fld}); shading interp; cbar=gcmfaces_cmap_cbar(cc/5); title(title0);
        %S
        subplot(2,1,2);
        [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.fldSzonmean(iLat,:,:)),'first');
        X=TT*ones(1,length(mygrid.RC)); Y=ones(length(TT),1)*(mygrid.RC'); X=X'; Y=Y';
        title0=['mean S minus first year (psu) at lat ~ ' num2str(listLatPlot(iLatPlot))];
        depthStretchPlot('pcolor',{X,Y,fld}); shading interp; cbar=gcmfaces_cmap_cbar(cc/25); title(title0);

        myCaption={'mean temperature (top; K) and salinity (bottom; psu) minus first year ',...
            ['at lat $\\approx$ ' num2str(listLatPlot(iLatPlot))]};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;


    if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'szmtime')));

    if addToTex; write2tex(fileTex,1,'zonal mean properties (surface)',2); end;

    %zonal mean SST/SSS:
    figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2]; kk=1;
    x=TT*ones(1,length(mygrid.LATS)); y=ones(nt,1)*mygrid.LATS'; x=x'; y=y';
    %T
    subplot(2,1,1);
    fld=squeeze(alldiag.fldTzonmean(:,kk,:)); [tmp1,fld]=annualmean(alldiag.listTimes,fld,'first');
    if ~doAnomalies; cc=cc0*3; else; cc=scaleAnom/10*cc0; end;
    pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]); cbar=gcmfaces_cmap_cbar(cc);
    title('zonal mean SST minus first year (degC)');
    %S
    subplot(2,1,2);
    fld=squeeze(alldiag.fldSzonmean(:,kk,:)); [tmp1,fld]=annualmean(alldiag.listTimes,fld,'first');
    if ~doAnomalies; cc=cc0*2/5; else; cc=scaleAnom/10*cc0/2; end;
    pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]); cbar=gcmfaces_cmap_cbar(cc);
    title('zonal mean SSS minus first year (psu)');
    %
    myCaption={'zonal mean temperature (degC; top) and salinity ',...
        '(psu; bottom) minus first year (psu)',[' at ' num2str(-mygrid.RC(kk)) 'm depth']};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2]; kk=1;
    if doAnomalies; cc=scaleAnom/10*cc/2; end;
    x=TT*ones(1,length(mygrid.LATS)); y=ones(nt,1)*mygrid.LATS'; x=x'; y=y';
    %SSH
    subplot(2,1,1);
    [tmp1,fld]=annualmean(alldiag.listTimes,alldiag.fldSSHzonmean,'first');
    pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
    cbar=gcmfaces_cmap_cbar(2*cc/25); title(['zonal mean SSH minus first year (INCLUDING ice, in m)']);
    %ETAN
    subplot(2,1,2);
    [tmp1,fld]=annualmean(alldiag.listTimes,alldiag.fldETANzonmean,'first');
    pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
    cbar=gcmfaces_cmap_cbar(4*cc/5); title(['zonal mean SSH minus first year (EXCLUDING ice, in m)']);
    %
    myCaption={'zonal mean SSH (m, uncorrected free surface) minus first year, including ice (top) and below ice (bottom) '};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2]; kk=1;
    if doAnomalies; cc=scaleAnom/10*cc; end;
    x=TT*ones(1,length(mygrid.LATS)); y=ones(nt,1)*mygrid.LATS'; x=x'; y=y';
    %SSH
    subplot(2,1,1);
    fld=squeeze(alldiag.fldSIzonmean(:,tt));
    pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
    if doAnomalies; cc2=scaleAnom/50*cc; else; cc2=0.5+cc/4; end;
    cbar=gcmfaces_cmap_cbar(cc2); title(['zonal mean Ice conc. -- in m ']);
    %MLD
    subplot(2,1,2);
    fld=squeeze(alldiag.fldMLDzonmean(:,tt));
    pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
    if doAnomalies; cc2=scaleAnom*40*cc; else; cc2=200+500*cc/5; end;
    cbar=gcmfaces_cmap_cbar(cc2); title(['zonal mean MLD -- in m ']);
    %
    myCaption={'zonal mean ice concentration (no units) and mixed layer depth (m)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    end;


    if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'icetime')));

    if addToTex; write2tex(fileTex,1,'seaice time series',2); end;

    %sea ice area
    figureL;
    subplot(2,1,1); vec=alldiag.IceAreaNorth/1e12;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 20]);
    if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
    grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice cover (in 10^{12}m^2)');
    subplot(2,1,2); vec=alldiag.IceAreaSouth/1e12;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 25]);
    grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice cover (in 10^{12}m^2)');
    if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
    myCaption={'sea ice cover (in $10^{12}m^2$) in northern (top) and southern (bottom) hemisphere'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %sea ice volume
    figureL;
    subplot(2,1,1); vec=alldiag.IceVolNorth/1e12;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 50]);
    grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice volume (in 10^{12m}^3)');
    if doAnomalies; aa(3:4)=scaleAnom*[-1 1]*2; axis(aa); end;
    subplot(2,1,2); vec=alldiag.IceVolSouth/1e12;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 15]);
    grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice volume (in 10^{12}m^3)');
    if doAnomalies; aa(3:4)=scaleAnom*[-1 1]*2; axis(aa); end;
    myCaption={'sea ice volume (in $10^{12}m^3$) in northern (top) and southern (bottom) hemisphere'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %snow volume
    figureL;
    subplot(2,1,1); vec=alldiag.SnowVolNorth/1e12;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 5]);
    grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere snow volume (in 10^{12}m^3)');
    if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
    subplot(2,1,2); vec=alldiag.SnowVolSouth/1e12;
    plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 4]);
    grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere snow volume (in 10^{12}m^3)');
    if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
    myCaption={'snow volume (in $10^{12}m^3$) in northern (top) and southern (bottom) hemisphere'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    if ~doAnomalies; %does not work for anomalies -> need to compute ratio in basic_diags_ecco
        %sea ice thickness
        figureL;
        subplot(2,1,1); vec=alldiag.IceVolNorth./alldiag.IceAreaNorth;
        plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 5]);
        grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice thickness (in m)');
        if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
        subplot(2,1,2); vec=alldiag.IceVolSouth./alldiag.IceAreaSouth;
        plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 2]);
        grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice thickness (in m)');
        if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
        myCaption={'sea ice thickness (in m) in northern (top) and southern (bottom) hemisphere'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        %sea ice thickness
        figureL;
        subplot(2,1,1); vec=alldiag.SnowVolNorth./alldiag.IceAreaNorth;
        plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 0.4]);
        grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere snow thickness (in m)');
        if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
        subplot(2,1,2); vec=alldiag.SnowVolSouth./alldiag.IceAreaSouth;
        plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 0.25]);
        grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere snow thickness (in m)');
        if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
        myCaption={'snow thickness (in m) in northern (top) and southern (bottom) hemisphere'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;

end;
1
2	if userStep==1;%diags to be computed
3	listDiags='fldTzonmean fldSzonmean fldETANzonmean fldSSHzonmean fldSIzonmean fldMLDzonmean';
4	listDiags=[listDiags ' IceAreaNorth IceAreaSouth IceVolNorth IceVolSouth SnowVolNorth SnowVolSouth'];
5	listDiags=[listDiags ' Ueq Teq SSHglo Tglo Sglo TgloProf SgloProf'];
6	elseif userStep==2;%input files and variables
7	listFlds={ 'THETA','SALT','ETAN','SIarea','SIheff','SIhsnow','sIceLoad','MXLDEPTH','UVELMASS','VVELMASS'};
8	listFldsNames=deblank(listFlds);
9	listFiles={'state_2d_set1','other_2d_set1','state_3d_set1','trsp_3d_set1'};
10	listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/'],[dirModel 'diags/']};
11	elseif userStep==3;%computational part;
12	%mask and re-arrange fields:
13	fldT=THETA.mygrid.mskC; fldS=SALT.mygrid.mskC;
14	fldETAN=ETAN.*mygrid.mskC(:,:,1);
15	fldSSH=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
16	fldSI=SIarea.*mygrid.mskC(:,:,1);
17	fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
18	ux=UVELMASS.mygrid.mskW; vy=VVELMASS.mygrid.mskS;
19	[ue,un]=calc_UEVNfromUXVY(ux,vy);
20
21	%compute zonal means:
22	[fldTzonmean]=calc_zonmean_T(fldT);
23	[fldSzonmean]=calc_zonmean_T(fldS);
24	[fldETANzonmean]=calc_zonmean_T(fldETAN);
25	[fldSSHzonmean]=calc_zonmean_T(fldSSH);
26	[fldSIzonmean]=calc_zonmean_T(fldSI);
27	[fldMLDzonmean]=calc_zonmean_T(fldMLD);
28
29	%compute hemispheric ice sums:
30	fld=SIarea.mygrid.RAC.(mygrid.YC>0); IceAreaNorth=nansum(fld);
31	fld=SIarea.mygrid.RAC.(mygrid.YC<0); IceAreaSouth=nansum(fld);
32	fld=SIheff.mygrid.RAC.(mygrid.YC>0); IceVolNorth=nansum(fld);
33	fld=SIheff.mygrid.RAC.(mygrid.YC<0); IceVolSouth=nansum(fld);
34	fld=SIhsnow.mygrid.RAC.(mygrid.YC>0); SnowVolNorth=nansum(fld);
35	fld=SIhsnow.mygrid.RAC.(mygrid.YC<0); SnowVolSouth=nansum(fld);
36
37	%equatorial sections of T and U:
38	[secX,secY,Teq]=gcmfaces_section([],0,fldT,1);
39	[secX,secY,Ueq]=gcmfaces_section([],0,ue,1);
40
41	%global means and profiles:
42	msk=mygrid.RAC.*mygrid.mskC(:,:,1);
43	SSHglo=nansum(fldSSH.*msk)./nansum(msk);
44	%
45	msk=mygrid.mskC.*mygrid.hFacC;
46	msk=msk.mk3D(mygrid.RAC,msk).mk3D(mygrid.DRF,msk);
47	Tglo=nansum(fldT.*msk)./nansum(msk);
48	Sglo=nansum(fldS.*msk)./nansum(msk);
49	nr=length(mygrid.RC);
50	TgloProf=zeros(nr,1); SgloProf=zeros(nr,1);
51	for kk=1:nr;
52	rac=mygrid.RAC.*mygrid.mskC(:,:,kk);
53	TgloProf(kk)=nansum(rac.*fldT(:,:,kk))/nansum(rac);
54	SgloProf(kk)=nansum(rac.*fldS(:,:,kk))/nansum(rac);
55	end;
56
57
58	%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
59	%===================== PLOTTING SEQUENCE BEGINS =========================%
60
61	elseif userStep==-1;%plotting
62
63	if isempty(setDiagsParams);
64	choicePlot={'all'};
65	else;
66	choicePlot=setDiagsParams;
67	end;
68
69	%ensure backward compatibility
70	if ~isfield(alldiag,'fldSSHzonmean'); alldiag.fldSSHzonmean=alldiag.fldETANLEADSzonmean; end;
71	if ~isfield(alldiag,'SSHglo'); alldiag.SSHglo=NaN*alldiag.Tglo; end;
72
73	%determine number of years in alldiag.listTimes
74	myTimes=alldiag.listTimes;
75	%determine the number of records in one year (lYear)
76	tmp1=mean(myTimes(2:end)-myTimes(1:end-1));
77	lYear=round(1/tmp1);
78	%in case when lYear<2 we use records as years
79	if ~(lYear>=2); lYear=1; myTimes=[1:length(myTimes)]; end;
80	%determine the number of full years (nYears)
81	nYears=floor(length(myTimes)/lYear);
82
83	if nYears>1&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'zmtend')));
84
85	if addToTex; write2tex(fileTex,1,'zonal mean tendencies',2); end;
86
87	figureL; set(gcf,'Renderer','zbuffer');
88	%cc=[-2:0.25:2];
89	cc=1/100*[[-200:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:200]];
90	if doAnomalies; cc=scaleAnom/10*cc; end;
91	X=mygrid.LATSones(1,length(mygrid.RC)); Y=ones(length(mygrid.LATS),1)(mygrid.RC');
92	%T
93	subplot(2,1,1);
94	fld=annualmean(alldiag.listTimes,alldiag.fldTzonmean,'last')-annualmean(alldiag.listTimes,alldiag.fldTzonmean,'first');
95	depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
96	cbar=gcmfaces_cmap_cbar(cc); title('zonal mean T anomaly');
97	%S
98	subplot(2,1,2);
99	fld=annualmean(alldiag.listTimes,alldiag.fldSzonmean,'last')-annualmean(alldiag.listTimes,alldiag.fldSzonmean,'first');
100	depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
101	cbar=gcmfaces_cmap_cbar(cc/4); title('zonal mean S anomaly');
102	%
103	myCaption={myYmeanTxt,', last year minus first year -- zonal mean temperature (degC; top) and salinity (psu; bottom)'};
104	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
105
106	end;
107
108	if sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'eq'));
109
110	if addToTex; write2tex(fileTex,1,'equatorial sections',2); end;
111
112	%time mean equator sections :
113	figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2];
114	[secX,secY,LONeq]=gcmfaces_section([],0,mygrid.XC,1);
115	X=LONeqones(1,length(mygrid.RC)); Y=ones(length(LONeq),1)mygrid.RC';
116	X=circshift(X,[-180 0]); X(X<0)=X(X<0)+360;
117	%T
118	subplot(2,1,1);
119	fld=annualmean(alldiag.listTimes,alldiag.Teq,'all'); fld=circshift(fld,[-180 0]);
120	depthStretchPlot('pcolor',{X,Y,fld},[0:50:350],[0 350 350]); shading interp;
121	if ~doAnomalies; cc=18+6cc0; else; cc=scaleAnom/10cc0*2; end;
122	cbar=gcmfaces_cmap_cbar(cc); title('equator T (degC)');
123	%U
124	subplot(2,1,2);
125	fld=annualmean(alldiag.listTimes,alldiag.Ueq,'all'); fld=circshift(fld,[-180 0]);
126	depthStretchPlot('pcolor',{X,Y,fld},[0:50:350],[0 350 350]); shading interp;
127	if ~doAnomalies; cc=cc0/52; else; cc=scaleAnom/10cc0/5; end;
128	cbar=gcmfaces_cmap_cbar(cc); title('equator zonal velocity (m/s)');
129	%
130	myCaption={myYmeanTxt,' mean -- equator temperature (degC;top) and zonal velocity (m/s;bottom)'};
131	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
132
133	end;
134
135	if multiTimes&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'gmtime')));
136
137	if addToTex; write2tex(fileTex,1,'global mean properties',2); end;
138
139	%global mean T/S:
140	figureL;
141	subplot(3,1,1); vec=alldiag.SSHglo;
142	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
143	aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
144	legend('monthly','annual mean'); title('Global Mean Sea level (m, uncorrected free surface)');
145	subplot(3,1,2); vec=alldiag.Tglo;
146	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
147	aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
148	legend('monthly','annual mean'); title('Global Mean Temperature (degC)');
149	subplot(3,1,3); vec=alldiag.Sglo;
150	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
151	aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
152	legend('monthly','annual mean'); title('Global Mean Salinity (psu)');
153	myCaption={'global mean T (degC; top) and S (psu; bottom)'};
154	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
155
156	%global mean T/S profiles:
157	figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2];
158	X=TTones(1,length(mygrid.RC)); Y=ones(length(TT),1)(mygrid.RC'); X=X'; Y=Y';
159	%T
160	subplot(2,1,1);
161	[tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.TgloProf),'first');
162	depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
163	if ~doAnomalies; cc=cc0/5; else; cc=scaleAnom/10*cc0/10; end;
164	cbar=gcmfaces_cmap_cbar(cc); title('global mean T minus first year (K)');
165	%S
166	subplot(2,1,2);
167	[tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.SgloProf),'first');
168	depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
169	if ~doAnomalies; cc=cc0/50; else; cc=scaleAnom/10*cc0/50; end;
170	cbar=gcmfaces_cmap_cbar(cc); title('global mean S minus first year (psu)');
171	%
172	myCaption={'global mean temperature (K; top) and salinity (psu; bottom) minus first year'};
173	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
174
175	end;
176
177
178	if multiTimes&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'lzmtime')));
179
180	if addToTex; write2tex(fileTex,1,'zonal mean properties',2); end;
181
182	listLatPlot=[-75 -65 -45 -25 0 25 45 65 75];
183	%zonal mean temperature profile at selected latitude:
184	for iLatPlot=1:length(listLatPlot);
185	tmp1=abs(mygrid.LATS-listLatPlot(iLatPlot));
186	iLat=find(tmp1==min(tmp1)); iLat=iLat(1);
187	figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2];
188	if doAnomalies; cc=scaleAnom/10*cc; end;
189	%T
190	subplot(2,1,1);
191	[tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.fldTzonmean(iLat,:,:)),'first');
192	X=TTones(1,length(mygrid.RC)); Y=ones(length(TT),1)(mygrid.RC'); X=X'; Y=Y';
193	title0=['mean T minus first year (K) at lat ~ ' num2str(listLatPlot(iLatPlot))];
194	depthStretchPlot('pcolor',{X,Y,fld}); shading interp; cbar=gcmfaces_cmap_cbar(cc/5); title(title0);
195	%S
196	subplot(2,1,2);
197	[tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.fldSzonmean(iLat,:,:)),'first');
198	X=TTones(1,length(mygrid.RC)); Y=ones(length(TT),1)(mygrid.RC'); X=X'; Y=Y';
199	title0=['mean S minus first year (psu) at lat ~ ' num2str(listLatPlot(iLatPlot))];
200	depthStretchPlot('pcolor',{X,Y,fld}); shading interp; cbar=gcmfaces_cmap_cbar(cc/25); title(title0);
201
202	myCaption={'mean temperature (top; K) and salinity (bottom; psu) minus first year ',...
203	['at lat $\\approx$ ' num2str(listLatPlot(iLatPlot))]};
204	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
205	end;
206
207	end;
208
209
210	if multiTimes&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'szmtime')));
211
212	if addToTex; write2tex(fileTex,1,'zonal mean properties (surface)',2); end;
213
214	%zonal mean SST/SSS:
215	figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2]; kk=1;
216	x=TTones(1,length(mygrid.LATS)); y=ones(nt,1)mygrid.LATS'; x=x'; y=y';
217	%T
218	subplot(2,1,1);
219	fld=squeeze(alldiag.fldTzonmean(:,kk,:)); [tmp1,fld]=annualmean(alldiag.listTimes,fld,'first');
220	if ~doAnomalies; cc=cc03; else; cc=scaleAnom/10cc0; end;
221	pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]); cbar=gcmfaces_cmap_cbar(cc);
222	title('zonal mean SST minus first year (degC)');
223	%S
224	subplot(2,1,2);
225	fld=squeeze(alldiag.fldSzonmean(:,kk,:)); [tmp1,fld]=annualmean(alldiag.listTimes,fld,'first');
226	if ~doAnomalies; cc=cc02/5; else; cc=scaleAnom/10cc0/2; end;
227	pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]); cbar=gcmfaces_cmap_cbar(cc);
228	title('zonal mean SSS minus first year (psu)');
229	%
230	myCaption={'zonal mean temperature (degC; top) and salinity ',...
231	'(psu; bottom) minus first year (psu)',[' at ' num2str(-mygrid.RC(kk)) 'm depth']};
232	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
233
234	figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2]; kk=1;
235	if doAnomalies; cc=scaleAnom/10*cc/2; end;
236	x=TTones(1,length(mygrid.LATS)); y=ones(nt,1)mygrid.LATS'; x=x'; y=y';
237	%SSH
238	subplot(2,1,1);
239	[tmp1,fld]=annualmean(alldiag.listTimes,alldiag.fldSSHzonmean,'first');
240	pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
241	cbar=gcmfaces_cmap_cbar(2*cc/25); title(['zonal mean SSH minus first year (INCLUDING ice, in m)']);
242	%ETAN
243	subplot(2,1,2);
244	[tmp1,fld]=annualmean(alldiag.listTimes,alldiag.fldETANzonmean,'first');
245	pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
246	cbar=gcmfaces_cmap_cbar(4*cc/5); title(['zonal mean SSH minus first year (EXCLUDING ice, in m)']);
247	%
248	myCaption={'zonal mean SSH (m, uncorrected free surface) minus first year, including ice (top) and below ice (bottom) '};
249	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
250
251	figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2]; kk=1;
252	if doAnomalies; cc=scaleAnom/10*cc; end;
253	x=TTones(1,length(mygrid.LATS)); y=ones(nt,1)mygrid.LATS'; x=x'; y=y';
254	%SSH
255	subplot(2,1,1);
256	fld=squeeze(alldiag.fldSIzonmean(:,tt));
257	pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
258	if doAnomalies; cc2=scaleAnom/50*cc; else; cc2=0.5+cc/4; end;
259	cbar=gcmfaces_cmap_cbar(cc2); title(['zonal mean Ice conc. -- in m ']);
260	%MLD
261	subplot(2,1,2);
262	fld=squeeze(alldiag.fldMLDzonmean(:,tt));
263	pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
264	if doAnomalies; cc2=scaleAnom40cc; else; cc2=200+500*cc/5; end;
265	cbar=gcmfaces_cmap_cbar(cc2); title(['zonal mean MLD -- in m ']);
266	%
267	myCaption={'zonal mean ice concentration (no units) and mixed layer depth (m)'};
268	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
269
270	end;
271
272
273	if multiTimes&(sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'icetime')));
274
275	if addToTex; write2tex(fileTex,1,'seaice time series',2); end;
276
277	%sea ice area
278	figureL;
279	subplot(2,1,1); vec=alldiag.IceAreaNorth/1e12;
280	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 20]);
281	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
282	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice cover (in 10^{12}m^2)');
283	subplot(2,1,2); vec=alldiag.IceAreaSouth/1e12;
284	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 25]);
285	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice cover (in 10^{12}m^2)');
286	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
287	myCaption={'sea ice cover (in $10^{12}m^2$) in northern (top) and southern (bottom) hemisphere'};
288	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
289
290	%sea ice volume
291	figureL;
292	subplot(2,1,1); vec=alldiag.IceVolNorth/1e12;
293	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 50]);
294	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice volume (in 10^{12m}^3)');
295	if doAnomalies; aa(3:4)=scaleAnom[-1 1]2; axis(aa); end;
296	subplot(2,1,2); vec=alldiag.IceVolSouth/1e12;
297	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 15]);
298	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice volume (in 10^{12}m^3)');
299	if doAnomalies; aa(3:4)=scaleAnom[-1 1]2; axis(aa); end;
300	myCaption={'sea ice volume (in $10^{12}m^3$) in northern (top) and southern (bottom) hemisphere'};
301	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
302
303	%snow volume
304	figureL;
305	subplot(2,1,1); vec=alldiag.SnowVolNorth/1e12;
306	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 5]);
307	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere snow volume (in 10^{12}m^3)');
308	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
309	subplot(2,1,2); vec=alldiag.SnowVolSouth/1e12;
310	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 4]);
311	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere snow volume (in 10^{12}m^3)');
312	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
313	myCaption={'snow volume (in $10^{12}m^3$) in northern (top) and southern (bottom) hemisphere'};
314	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
315
316	if ~doAnomalies; %does not work for anomalies -> need to compute ratio in basic_diags_ecco
317	%sea ice thickness
318	figureL;
319	subplot(2,1,1); vec=alldiag.IceVolNorth./alldiag.IceAreaNorth;
320	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 5]);
321	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice thickness (in m)');
322	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
323	subplot(2,1,2); vec=alldiag.IceVolSouth./alldiag.IceAreaSouth;
324	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 2]);
325	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice thickness (in m)');
326	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
327	myCaption={'sea ice thickness (in m) in northern (top) and southern (bottom) hemisphere'};
328	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
329
330	%sea ice thickness
331	figureL;
332	subplot(2,1,1); vec=alldiag.SnowVolNorth./alldiag.IceAreaNorth;
333	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 0.4]);
334	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere snow thickness (in m)');
335	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
336	subplot(2,1,2); vec=alldiag.SnowVolSouth./alldiag.IceAreaSouth;
337	plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 0.25]);
338	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere snow thickness (in m)');
339	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
340	myCaption={'snow thickness (in m) in northern (top) and southern (bottom) hemisphere'};
341	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
342	end;
343
344	end;
345
346	end;