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	gforget	1.1
2			if userStep==1;%diags to be computed
3	gforget	1.4	listDiags='fldTzonmean fldSzonmean fldETANzonmean fldSSHzonmean fldSIzonmean fldMLDzonmean';
4	gforget	1.1	listDiags=[listDiags ' IceAreaNorth IceAreaSouth IceVolNorth IceVolSouth SnowVolNorth SnowVolSouth'];
5	gforget	1.4	listDiags=[listDiags ' Ueq Teq SSHglo Tglo Sglo TgloProf SgloProf'];
6	gforget	1.1	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	gforget	1.9	listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/'],[dirModel 'diags/']};
11	gforget	1.1	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	gforget	1.4	fldSSH=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
16	gforget	1.1	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	gforget	1.4	[fldSSHzonmean]=calc_zonmean_T(fldSSH);
26	gforget	1.1	[fldSIzonmean]=calc_zonmean_T(fldSI);
27			[fldMLDzonmean]=calc_zonmean_T(fldMLD);
28	gforget	1.10
29	gforget	1.1	%compute hemispheric ice sums:
30	gforget	1.10	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	gforget	1.1
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	gforget	1.4	msk=mygrid.RAC.*mygrid.mskC(:,:,1);
43			SSHglo=nansum(fldSSH.*msk)./nansum(msk);
44			%
45	gforget	1.1	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	gforget	1.2
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	gforget	1.4	%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	gforget	1.10	%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	gforget	1.2
85	gforget	1.4	if addToTex; write2tex(fileTex,1,'zonal mean tendencies',2); end;
86
87	gforget	1.2	figureL; set(gcf,'Renderer','zbuffer');
88			%cc=[-2:0.25:2];
89	gforget	1.8	cc=1/100*[[-200:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:200]];
90	gforget	1.2	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	gforget	1.4	if addToTex; write2tex(fileTex,1,'equatorial sections',2); end;
111
112	gforget	1.2	%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	gforget	1.4	if addToTex; write2tex(fileTex,1,'global mean properties',2); end;
138
139	gforget	1.2	%global mean T/S:
140			figureL;
141	gforget	1.4	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	gforget	1.7	legend('monthly','annual mean'); title('Global Mean Sea level (m, uncorrected free surface)');
145	gforget	1.5	subplot(3,1,2); vec=alldiag.Tglo;
146	gforget	1.2	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	gforget	1.5	subplot(3,1,3); vec=alldiag.Sglo;
150	gforget	1.2	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	gforget	1.4	if addToTex; write2tex(fileTex,1,'zonal mean properties',2); end;
181
182	gforget	1.8	listLatPlot=[-75 -65 -45 -25 0 25 45 65 75];
183	gforget	1.2	%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	gforget	1.4	if addToTex; write2tex(fileTex,1,'zonal mean properties (surface)',2); end;
213
214	gforget	1.2	%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	gforget	1.4	[tmp1,fld]=annualmean(alldiag.listTimes,alldiag.fldSSHzonmean,'first');
240	gforget	1.2	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	gforget	1.7	myCaption={'zonal mean SSH (m, uncorrected free surface) minus first year, including ice (top) and below ice (bottom) '};
249	gforget	1.2	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	gforget	1.3	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	gforget	1.2	%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	gforget	1.3	if doAnomalies; cc2=scaleAnom40cc; else; cc2=200+500*cc/5; end;
265			cbar=gcmfaces_cmap_cbar(cc2); title(['zonal mean MLD -- in m ']);
266	gforget	1.2	%
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	gforget	1.4	if addToTex; write2tex(fileTex,1,'seaice time series',2); end;
276
277	gforget	1.2	%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	heimbach	1.6	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice cover (in 10^{12}m^2)');
283	gforget	1.2	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	heimbach	1.6	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice cover (in 10^{12}m^2)');
286	gforget	1.2	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
287	heimbach	1.6	myCaption={'sea ice cover (in $10^{12}m^2$) in northern (top) and southern (bottom) hemisphere'};
288	gforget	1.2	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	heimbach	1.6	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice volume (in 10^{12m}^3)');
295	gforget	1.2	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	heimbach	1.6	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice volume (in 10^{12}m^3)');
299	gforget	1.2	if doAnomalies; aa(3:4)=scaleAnom[-1 1]2; axis(aa); end;
300	heimbach	1.6	myCaption={'sea ice volume (in $10^{12}m^3$) in northern (top) and southern (bottom) hemisphere'};
301	gforget	1.2	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	heimbach	1.6	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere snow volume (in 10^{12}m^3)');
308	gforget	1.2	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	heimbach	1.6	grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere snow volume (in 10^{12}m^3)');
312	gforget	1.2	if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
313	heimbach	1.6	myCaption={'snow volume (in $10^{12}m^3$) in northern (top) and southern (bottom) hemisphere'};
314	gforget	1.2	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	gforget	1.1	end;