matlab_class/gcmfaces_diags/diags_set_B.m


doOmit3Dfields=0;
if mygrid.memoryLimit~=0; doOmit3Dfields=1; end;

if userStep==1;%diags to be computed
    %here _s stands for cumulative sum and _2s for cumulative sum of squares
    %the _s should be stated first
    listDiags=['SIatmQnt_s SIatmFW_s oceQnet_s oceFWflx_s ' ...
               'fldTZ_s fldTM_s curlTau_s fldETAN_s fldETANLEADS_s fldMLD_s ' ...
               'SIatmQnt_2s SIatmFW_2s oceQnet_2s oceFWflx_2s ' ...
               'fldTZ_2s fldTM_2s curlTau_2s fldETAN_2s fldETANLEADS_2s fldMLD_2s'];
    if ~doOmit3Dfields;
      listDiags=[listDiags 'THETA_s SALT_s WVELMASS_s THETA_2s SALT_2s WVELMASS_2s'];
    end;

elseif userStep==2;%input files and variables
    listFlds={  'SIatmQnt','SIatmFW ','oceQnet ','oceFWflx','oceTAUX','oceTAUY','ETAN','sIceLoad','MXLDEPTH'};
    if ~doOmit3Dfields;
      listFlds={listFlds{:} ,'THETA','SALT','WVELMASS'};
    end;
    listFldsNames=deblank(listFlds);
    listFiles={'state_3d_set1','trsp_3d_set1','state_2d_set1','other_2d_set1'};
    listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/' ]};

elseif userStep==3;%computation
    %mask fields:
    SIatmQnt=SIatmQnt.*mygrid.mskC(:,:,1);
    SIatmFW=SIatmFW.*mygrid.mskC(:,:,1);
    oceQnet=oceQnet.*mygrid.mskC(:,:,1);
    oceFWflx=oceFWflx.*mygrid.mskC(:,:,1);
    fldTX=oceTAUX.*mygrid.mskW(:,:,1);
    fldTY=oceTAUY.*mygrid.mskS(:,:,1);
    %compute Eastward/Northward wind stresses:
    [fldTZ,fldTM]=calc_UEVNfromUXVY(fldTX,fldTY);
    %compute wind stress curl:
    %curlTau=calc_UV_curl(fldTX, fldTY,1 );%the doMask argument should not matter as msk was already applied
    curlTau=NaN*fldTZ;
    %mask and re-arrange fields:
    fldETAN=ETAN.*mygrid.mskC(:,:,1);
    fldETANLEADS=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
    fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
    %
    if ~doOmit3Dfields;
      THETA=THETA.*mygrid.mskC;
      SALT=SALT.*mygrid.mskC;
      WVELMASS=WVELMASS.*mygrid.mskC;
    end;
    %
    if ii>1;
      fileMatPrev=['diags_set_' tmp1 '_' num2str(listTimes(ii-1)) '.mat'];
      listTimesBak=listTimes;
      load([dirMat fileMatPrev]);
      listTimes=listTimesBak;
    end;
    for jj=1:length(listDiags)/2;
        myDiag=listDiags{jj}(1:end-2);
        if ii==1;
            eval([myDiag '_s=0*' myDiag ';']);
            eval([myDiag '_2s=0*' myDiag '.^2;']);
        end;
        eval([myDiag '_s=' myDiag '_s+' myDiag ';']);
        eval([myDiag '_2s=' myDiag '_2s+' myDiag '.^2;']);
    end;

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

elseif userStep==0;%loading / post-processing of mat files

  %load last cumsum
  alldiag=alldiag_load(dirMat,[fileMat '_*.mat'],'',myparms.recInAve(2));
  %load first cumsum
  if myparms.recInAve(1)>1;
    tmpdiag=alldiag_load(dirMat,[fileMat '_*.mat'],'',myparms.recInAve(1)-1);
    for ii=1:length(alldiag.listDiags);
         tmp0=alldiag.listDiags{ii};
         if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps');
           tmp1=getfield(alldiag,alldiag.listDiags{ii});
           tmp2=getfield(tmpdiag,alldiag.listDiags{ii});
           alldiag=setfield(alldiag,tmp0,tmp1-tmp2);
         end;
    end;
  end;
  %accomodate missing fields (old MITgcm version)
  listFlds={'SIatmQnt_s','oceQnet_s','SIatmFW_s','oceFWflx_s',...
            'THETA_s','SALT_s','WVELMASS_s',...
            'fldTZ_s','fldTM_s','curlTau_s','fldETAN_s','fldETANLEADS_s',...
            'SIatmQnt_2s','oceQnet_2s','SIatmFW_2s','oceFWflx_2s',...
            'fldTZ_2s','fldTM_2s','curlTau_2s','fldETAN_2s','fldETANLEADS_2s',...
            'THETA_2s','SALT_2s','WVELMASS_2s'};
  for ii=1:length(listFlds);
    if ~sum(strcmp(fieldnames(alldiag),listFlds{ii}));
      eval(['alldiag.' listFlds{ii} '=NaN*alldiag.fldETAN_s;']);
    end;
  end;
  %
  n=diff(myparms.recInAve)+1;
  for ii=1:length(alldiag.listDiags)/2;
     tmp0=alldiag.listDiags{ii};
     if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps');
       tmp1=1/n*getfield(alldiag,tmp0);
       tmp2=1/n*getfield(alldiag,[tmp0(1:end-1) '2s']);
       tmp2=(tmp2-tmp1.^2);
       tmp2=n/(n-1)*tmp2;
       %tmp2(tmp2<0)=0;
       tmp2=sqrt(tmp2);
       alldiag=setfield(alldiag,[tmp0(1:end-2) '_mean'],tmp1);
       alldiag=setfield(alldiag,[tmp0(1:end-2) '_std'],tmp2);
     end;
  end;

  diagsWereLoaded=1

elseif userStep==-1;%plotting

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


%===== start with state variables

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

    if addToTex; write2tex(fileTex,1,'sea surface height',2); end;

    %ETAN:
    fld=alldiag.fldETAN_mean;
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/100; title0='sea surface height (EXCLUDING ice)';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.05; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- sea surface height (EXCLUDING ice, in m)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %ETANLEADS:
    fld=alldiag.fldETANLEADS_mean;
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/100; title0='sea surface height (INCLUDING ice)';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.05; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- sea surface height (INCLUDING ice, in m)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    if multiTimes;
        %ETAN:
        fld=alldiag.fldETAN_std;
        cc=[0:25:500]/2500; title0='std(ETAN)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.02; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,' standard deviation -- sea surface height (EXCLUDING ice, in m)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        %ETANLEADS:
        fld=alldiag.fldETANLEADS_std;
        cc=[0:25:500]/2500; title0='std(ETANLEADS)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.02; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,' standard deviation -- sea surface height (INCLUDING ice, in m)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;

  if ( sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'subsrfc')) )&(~doOmit3Dfields);

    if addToTex; write2tex(fileTex,1,'3D state variables',2); end;

    cc_mean=1/100*[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]];
    cc_std=1/250*[0:25:500];
    kkList=[1 11 20 28 37 44];
    vvList={'THETA','SALT','WVELMASS'};

    for vv=1:length(vvList);
      vvtxt=vvList{vv};
      eval(['fld_mean=alldiag.' vvtxt '_mean;']);
      eval(['fld_std=alldiag.' vvtxt '_std;']);
      if strcmp(vvtxt,'WVELMASS');%convert into mm/day
           fld_mean=fld_mean*86400*365/1e3;
           fld_std=fld_std*86400*365/1e3;
      end;
      for kk=kkList;
        if strcmp(vvtxt,'WVELMASS')&kk==1; kk=2; end;
        if mygrid.RC(kk)>=-50; facD=0.5;
        elseif mygrid.RC(kk)>=-150; facD=1;
        elseif mygrid.RC(kk)>=-500; facD=2.5;
        elseif mygrid.RC(kk)>=-2000; facD=10;
        else; facD=20;
        end;
        %
        if strcmp(vvtxt,'THETA'); nm='temperature (in degrees)'; facV=1;
        elseif strcmp(vvtxt,'SALT'); nm='salinity (in psu)'; facV=10;
        elseif strcmp(vvtxt,'WVELMASS'); nm='vertical velocity (in mm/year)'; facV=10;
        else; error('not yet implemented');
        end;
        %

    title0=sprintf('%s at %4dm',nm,round(-mygrid.RC(kk)));

    %time mean map:
    fld=fld_mean(:,:,kk).*mygrid.mskC(:,:,kk);
    if doAnomalies;
      cc=scaleAnom/facV/facD*cc_mean; 
    elseif strcmp(vvtxt,'WVELMASS');
      cc=1/facV*cc_mean; 
    else;
      cc=prctile(facV*fld,[2.5 97.5]);
      cc=[floor(cc(1)):ceil(cc(end))]/facV;
      %ensure at least 10 color levels
      fac0=1;
      while length(cc)<15;
        fac0=fac0*2;
        cc=prctile(fac0*facV*fld,[2.5 97.5]);
        cc=[floor(cc(1)):ceil(cc(end))]/facV/fac0;
      end;
      %ensure at most 10 color levels
      fac0=1;
      while length(cc)>30;
        fac0=fac0/2;
        cc=prctile(fac0*facV*fld,[2.5 97.5]);
        cc=[floor(cc(1)):ceil(cc(end))]/facV/fac0;
      end;
    end;
    %
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- ',title0};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    if 0;%multiTimes;
        fld=fld_std(:,:,kk).*mygrid.mskC(:,:,kk);
        cc=1/facV/facD*cc_std; 
        if strcmp(vvtxt,'WVELMASS'); cc=1/facV*cc_std; end;
        if doAnomalies; cc=scaleAnom*cc; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,' standard deviation -- ',title0};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;%for kk=...
    end;%for vv=...
    end;%if ( sum(strcmp(...


%now do surface fluxes and forcing fields

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

    if addToTex; write2tex(fileTex,1,'air-sea heat flux',2); end;

    %qnet from ocean+ice:
    fld=-alldiag.SIatmQnt_mean;
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]; title0='QNET to ocean+ice';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*10; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- QNET to ocean+ice (W/m2)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %qnet to ocean:
    fld=alldiag.oceQnet_mean;
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]; title0='QNET to ocean';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*10; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- QNET to ocean (W/m2)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
   
    if multiTimes;
        %qnet from ocean+ice:
        fld=alldiag.SIatmQnt_std;
        cc=[[0:5:25] 35 [50:25:200] [250 300]]; title0='std(QNET to ocean+ice)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*5; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,'  standard deviation -- QNET to ocean+ice (W/m2)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        %qnet from ocean:
        fld=alldiag.oceQnet_std;
        cc=[[0:5:25] 35 [50:25:200] [250 300]]; title0='std(QNET to ocean)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*5; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,'  standard deviation -- QNET to ocean (W/m2)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;


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

    if addToTex; write2tex(fileTex,1,'air-sea fresh water flux',2); end;

    %FW flux from ocean+ice:
    fld=-alldiag.SIatmFW_mean/1000;%conversion to m/s
    fld=fld*86400*1000;%conversion to mm/day
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]*0.06; title0='EMPMR to ocean+ice';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.5; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- EMPMR from ocean+ice (mm/day)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %FW flux from ocean:
    fld=-alldiag.oceFWflx_mean/1000;%conversion to m/s
    fld=fld*86400*1000;%conversion to mm/day
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]*0.06; title0='EMPMR to ocean';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.5; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- EMPMR from ocean (mm/day)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    if multiTimes;
        %empmr from ocean+ice:
        fld=alldiag.SIatmFW_std*86400;%conversion to mm/day
        cc=[[0:5:25] 35 [50:25:200] [250 300]]*0.04; title0='std(EMPMR to ocean+ice)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.5; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,' standard deviation -- EMPMR to ocean+ice (W/m2)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        %empmr from ocean:
        fld=alldiag.oceFWflx_std*86400;%conversion to mm/day
        cc=[[0:5:25] 35 [50:25:200] [250 300]]*0.04; title0='std(EMPMR to ocean)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.5; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,' standard deviation -- EMPMR to ocean (W/m2)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;

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

    if addToTex; write2tex(fileTex,1,'surface wind stress',2); end;

    %zonal wind stress:
    fld=alldiag.fldTZ_mean;
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='zonal wind stress';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- zonal wind stress (N/m2)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %meridional wind stress:
    fld=alldiag.fldTM_mean;
    cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='meridional wind stress';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
    figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    myCaption={myYmeanTxt,'mean -- meridional wind stress (N/m2)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %fld=alldiag.curlTau_mean;
    %cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/5e8; title0='wind stress curl';
    %if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
    %figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
    %myCaption={myYmeanTxt,'mean -- wind stress curl (N/m3)'};
    %if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    if multiTimes;
        %zonal wind stress:
        fld=alldiag.fldTZ_std;
        cc=[[0:5:25] 35 [50:25:200] [250 300]]/2000; title0='std(tauZ)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.005; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,'  standard deviation -- tauZ (W/m2)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        %meridional wind stress:
        fld=alldiag.fldTM_std;
        cc=[[0:5:25] 35 [50:25:200] [250 300]]/2000; title0='std(tauM)';
        if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.005; end;
        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        myCaption={myYmeanTxt,' standard deviation -- tauM (W/m2)'};
        if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

        %fld=alldiag.curlTau_std;
        %cc=[[0:5:25] 35 [50:25:200] [250 300]]/1e9; title0='wind stress curl';
        %if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
        %figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
        %myCaption={myYmeanTxt,'standard deviation -- tauCurl (N/m3)'};
        %if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;

end;


1	gforget	1.1
2	gforget	1.5	doOmit3Dfields=0;
3	gforget	1.6	if mygrid.memoryLimit~=0; doOmit3Dfields=1; end;
4	gforget	1.5
5	gforget	1.1	if userStep==1;%diags to be computed
6	gforget	1.2	%here _s stands for cumulative sum and _2s for cumulative sum of squares
7			%the _s should be stated first
8	gforget	1.5	listDiags=['SIatmQnt_s SIatmFW_s oceQnet_s oceFWflx_s ' ...
9	gforget	1.2	'fldTZ_s fldTM_s curlTau_s fldETAN_s fldETANLEADS_s fldMLD_s ' ...
10	gforget	1.5	'SIatmQnt_2s SIatmFW_2s oceQnet_2s oceFWflx_2s ' ...
11	gforget	1.2	'fldTZ_2s fldTM_2s curlTau_2s fldETAN_2s fldETANLEADS_2s fldMLD_2s'];
12	gforget	1.5	if ~doOmit3Dfields;
13			listDiags=[listDiags 'THETA_s SALT_s WVELMASS_s THETA_2s SALT_2s WVELMASS_2s'];
14			end;
15	gforget	1.4
16	gforget	1.1	elseif userStep==2;%input files and variables
17	gforget	1.5	listFlds={ 'SIatmQnt','SIatmFW ','oceQnet ','oceFWflx','oceTAUX','oceTAUY','ETAN','sIceLoad','MXLDEPTH'};
18			if ~doOmit3Dfields;
19			listFlds={listFlds{:} ,'THETA','SALT','WVELMASS'};
20			end;
21	gforget	1.1	listFldsNames=deblank(listFlds);
22	gforget	1.5	listFiles={'state_3d_set1','trsp_3d_set1','state_2d_set1','other_2d_set1'};
23			listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/' ]};
24	gforget	1.4
25			elseif userStep==3;%computation
26	gforget	1.1	%mask fields:
27			SIatmQnt=SIatmQnt.*mygrid.mskC(:,:,1);
28			SIatmFW=SIatmFW.*mygrid.mskC(:,:,1);
29			oceQnet=oceQnet.*mygrid.mskC(:,:,1);
30			oceFWflx=oceFWflx.*mygrid.mskC(:,:,1);
31			fldTX=oceTAUX.*mygrid.mskW(:,:,1);
32			fldTY=oceTAUY.*mygrid.mskS(:,:,1);
33			%compute Eastward/Northward wind stresses:
34			[fldTZ,fldTM]=calc_UEVNfromUXVY(fldTX,fldTY);
35			%compute wind stress curl:
36	gforget	1.6	%curlTau=calc_UV_curl(fldTX, fldTY,1 );%the doMask argument should not matter as msk was already applied
37			curlTau=NaN*fldTZ;
38	gforget	1.1	%mask and re-arrange fields:
39			fldETAN=ETAN.*mygrid.mskC(:,:,1);
40			fldETANLEADS=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
41			fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
42			%
43	gforget	1.5	if ~doOmit3Dfields;
44			THETA=THETA.*mygrid.mskC;
45			SALT=SALT.*mygrid.mskC;
46			WVELMASS=WVELMASS.*mygrid.mskC;
47			end;
48	gforget	1.2	%
49			if ii>1;
50			fileMatPrev=['diags_set_' tmp1 '_' num2str(listTimes(ii-1)) '.mat'];
51	gforget	1.4	listTimesBak=listTimes;
52	gforget	1.2	load([dirMat fileMatPrev]);
53	gforget	1.4	listTimes=listTimesBak;
54	gforget	1.2	end;
55			for jj=1:length(listDiags)/2;
56			myDiag=listDiags{jj}(1:end-2);
57			if ii==1;
58			eval([myDiag '_s=0*' myDiag ';']);
59			eval([myDiag '_2s=0*' myDiag '.^2;']);
60			end;
61			eval([myDiag '_s=' myDiag '_s+' myDiag ';']);
62			eval([myDiag '_2s=' myDiag '_2s+' myDiag '.^2;']);
63			end;
64
65	gforget	1.4	%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
66			%===================== PLOTTING SEQUENCE BEGINS =========================%
67
68			elseif userStep==0;%loading / post-processing of mat files
69
70			%load last cumsum
71			alldiag=alldiag_load(dirMat,[fileMat '_*.mat'],'',myparms.recInAve(2));
72			%load first cumsum
73			if myparms.recInAve(1)>1;
74			tmpdiag=alldiag_load(dirMat,[fileMat '_*.mat'],'',myparms.recInAve(1)-1);
75			for ii=1:length(alldiag.listDiags);
76			tmp0=alldiag.listDiags{ii};
77			if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps');
78			tmp1=getfield(alldiag,alldiag.listDiags{ii});
79			tmp2=getfield(tmpdiag,alldiag.listDiags{ii});
80			alldiag=setfield(alldiag,tmp0,tmp1-tmp2);
81			end;
82			end;
83			end;
84			%accomodate missing fields (old MITgcm version)
85			listFlds={'SIatmQnt_s','oceQnet_s','SIatmFW_s','oceFWflx_s',...
86	gforget	1.5	'THETA_s','SALT_s','WVELMASS_s',...
87	gforget	1.4	'fldTZ_s','fldTM_s','curlTau_s','fldETAN_s','fldETANLEADS_s',...
88			'SIatmQnt_2s','oceQnet_2s','SIatmFW_2s','oceFWflx_2s',...
89	gforget	1.5	'fldTZ_2s','fldTM_2s','curlTau_2s','fldETAN_2s','fldETANLEADS_2s',...
90			'THETA_2s','SALT_2s','WVELMASS_2s'};
91	gforget	1.4	for ii=1:length(listFlds);
92			if ~sum(strcmp(fieldnames(alldiag),listFlds{ii}));
93			eval(['alldiag.' listFlds{ii} '=NaN*alldiag.fldETAN_s;']);
94			end;
95			end;
96			%
97			n=diff(myparms.recInAve)+1;
98			for ii=1:length(alldiag.listDiags)/2;
99			tmp0=alldiag.listDiags{ii};
100			if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps');
101			tmp1=1/n*getfield(alldiag,tmp0);
102			tmp2=1/n*getfield(alldiag,[tmp0(1:end-1) '2s']);
103			tmp2=(tmp2-tmp1.^2);
104			tmp2=n/(n-1)*tmp2;
105			%tmp2(tmp2<0)=0;
106			tmp2=sqrt(tmp2);
107			alldiag=setfield(alldiag,[tmp0(1:end-2) '_mean'],tmp1);
108			alldiag=setfield(alldiag,[tmp0(1:end-2) '_std'],tmp2);
109			end;
110			end;
111
112			diagsWereLoaded=1
113
114			elseif userStep==-1;%plotting
115
116			if isempty(setDiagsParams);
117			choicePlot={'all'};
118			else;
119			choicePlot=setDiagsParams;
120			end;
121
122	gforget	1.6
123			%===== start with state variables
124
125			if sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'surface'));
126
127			if addToTex; write2tex(fileTex,1,'sea surface height',2); end;
128
129			%ETAN:
130			fld=alldiag.fldETAN_mean;
131			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/100; title0='sea surface height (EXCLUDING ice)';
132			if doAnomalies; cc=scaleAnom[-1:0.1:1]0.05; end;
133			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
134			myCaption={myYmeanTxt,'mean -- sea surface height (EXCLUDING ice, in m)'};
135			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
136
137			%ETANLEADS:
138			fld=alldiag.fldETANLEADS_mean;
139			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/100; title0='sea surface height (INCLUDING ice)';
140			if doAnomalies; cc=scaleAnom[-1:0.1:1]0.05; end;
141			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
142			myCaption={myYmeanTxt,'mean -- sea surface height (INCLUDING ice, in m)'};
143			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
144
145			if multiTimes;
146			%ETAN:
147			fld=alldiag.fldETAN_std;
148			cc=[0:25:500]/2500; title0='std(ETAN)';
149			if doAnomalies; cc=scaleAnom[0:0.1:1]0.02; end;
150			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
151			myCaption={myYmeanTxt,' standard deviation -- sea surface height (EXCLUDING ice, in m)'};
152			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
153
154			%ETANLEADS:
155			fld=alldiag.fldETANLEADS_std;
156			cc=[0:25:500]/2500; title0='std(ETANLEADS)';
157			if doAnomalies; cc=scaleAnom[0:0.1:1]0.02; end;
158			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
159			myCaption={myYmeanTxt,' standard deviation -- sea surface height (INCLUDING ice, in m)'};
160			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
161			end;
162
163			end;
164
165			if ( sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'subsrfc')) )&(~doOmit3Dfields);
166
167			if addToTex; write2tex(fileTex,1,'3D state variables',2); end;
168
169			cc_mean=1/100*[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]];
170			cc_std=1/250*[0:25:500];
171			kkList=[1 11 20 28 37 44];
172			vvList={'THETA','SALT','WVELMASS'};
173
174			for vv=1:length(vvList);
175			vvtxt=vvList{vv};
176			eval(['fld_mean=alldiag.' vvtxt '_mean;']);
177			eval(['fld_std=alldiag.' vvtxt '_std;']);
178			if strcmp(vvtxt,'WVELMASS');%convert into mm/day
179			fld_mean=fld_mean86400365/1e3;
180			fld_std=fld_std86400365/1e3;
181			end;
182			for kk=kkList;
183			if strcmp(vvtxt,'WVELMASS')&kk==1; kk=2; end;
184			if mygrid.RC(kk)>=-50; facD=0.5;
185			elseif mygrid.RC(kk)>=-150; facD=1;
186			elseif mygrid.RC(kk)>=-500; facD=2.5;
187			elseif mygrid.RC(kk)>=-2000; facD=10;
188			else; facD=20;
189			end;
190			%
191			if strcmp(vvtxt,'THETA'); nm='temperature (in degrees)'; facV=1;
192			elseif strcmp(vvtxt,'SALT'); nm='salinity (in psu)'; facV=10;
193			elseif strcmp(vvtxt,'WVELMASS'); nm='vertical velocity (in mm/year)'; facV=10;
194			else; error('not yet implemented');
195			end;
196			%
197
198			title0=sprintf('%s at %4dm',nm,round(-mygrid.RC(kk)));
199
200			%time mean map:
201			fld=fld_mean(:,:,kk).*mygrid.mskC(:,:,kk);
202			if doAnomalies;
203			cc=scaleAnom/facV/facD*cc_mean;
204			elseif strcmp(vvtxt,'WVELMASS');
205			cc=1/facV*cc_mean;
206			else;
207			cc=prctile(facV*fld,[2.5 97.5]);
208			cc=[floor(cc(1)):ceil(cc(end))]/facV;
209			%ensure at least 10 color levels
210			fac0=1;
211			while length(cc)<15;
212			fac0=fac0*2;
213			cc=prctile(fac0facVfld,[2.5 97.5]);
214			cc=[floor(cc(1)):ceil(cc(end))]/facV/fac0;
215			end;
216			%ensure at most 10 color levels
217			fac0=1;
218			while length(cc)>30;
219			fac0=fac0/2;
220			cc=prctile(fac0facVfld,[2.5 97.5]);
221			cc=[floor(cc(1)):ceil(cc(end))]/facV/fac0;
222			end;
223			end;
224			%
225			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
226			myCaption={myYmeanTxt,'mean -- ',title0};
227			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
228
229			if 0;%multiTimes;
230			fld=fld_std(:,:,kk).*mygrid.mskC(:,:,kk);
231			cc=1/facV/facD*cc_std;
232			if strcmp(vvtxt,'WVELMASS'); cc=1/facV*cc_std; end;
233			if doAnomalies; cc=scaleAnom*cc; end;
234			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
235			myCaption={myYmeanTxt,' standard deviation -- ',title0};
236			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
237			end;
238
239			end;%for kk=...
240			end;%for vv=...
241			end;%if ( sum(strcmp(...
242
243
244			%now do surface fluxes and forcing fields
245
246	gforget	1.4	if sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'qnet'));
247
248	gforget	1.6	if addToTex; write2tex(fileTex,1,'air-sea heat flux',2); end;
249
250	gforget	1.4	%qnet from ocean+ice:
251			fld=-alldiag.SIatmQnt_mean;
252			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]; title0='QNET to ocean+ice';
253			if doAnomalies; cc=scaleAnom[-1:0.1:1]10; end;
254			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
255			myCaption={myYmeanTxt,'mean -- QNET to ocean+ice (W/m2)'};
256			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
257
258			%qnet to ocean:
259			fld=alldiag.oceQnet_mean;
260			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]; title0='QNET to ocean';
261			if doAnomalies; cc=scaleAnom[-1:0.1:1]10; end;
262			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
263			myCaption={myYmeanTxt,'mean -- QNET to ocean (W/m2)'};
264			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
265
266			if multiTimes;
267			%qnet from ocean+ice:
268			fld=alldiag.SIatmQnt_std;
269			cc=[[0:5:25] 35 [50:25:200] [250 300]]; title0='std(QNET to ocean+ice)';
270			if doAnomalies; cc=scaleAnom[0:0.1:1]5; end;
271			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
272			myCaption={myYmeanTxt,' standard deviation -- QNET to ocean+ice (W/m2)'};
273			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
274
275			%qnet from ocean:
276			fld=alldiag.oceQnet_std;
277			cc=[[0:5:25] 35 [50:25:200] [250 300]]; title0='std(QNET to ocean)';
278			if doAnomalies; cc=scaleAnom[0:0.1:1]5; end;
279			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
280			myCaption={myYmeanTxt,' standard deviation -- QNET to ocean (W/m2)'};
281			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
282			end;
283
284			end;
285
286
287			if sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'fwf'));
288
289	gforget	1.6	if addToTex; write2tex(fileTex,1,'air-sea fresh water flux',2); end;
290
291	gforget	1.4	%FW flux from ocean+ice:
292			fld=-alldiag.SIatmFW_mean/1000;%conversion to m/s
293			fld=fld864001000;%conversion to mm/day
294			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]*0.06; title0='EMPMR to ocean+ice';
295			if doAnomalies; cc=scaleAnom[-1:0.1:1]0.5; end;
296			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
297			myCaption={myYmeanTxt,'mean -- EMPMR from ocean+ice (mm/day)'};
298			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
299
300			%FW flux from ocean:
301			fld=-alldiag.oceFWflx_mean/1000;%conversion to m/s
302			fld=fld864001000;%conversion to mm/day
303			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]*0.06; title0='EMPMR to ocean';
304			if doAnomalies; cc=scaleAnom[-1:0.1:1]0.5; end;
305			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
306			myCaption={myYmeanTxt,'mean -- EMPMR from ocean (mm/day)'};
307			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
308
309			if multiTimes;
310			%empmr from ocean+ice:
311			fld=alldiag.SIatmFW_std*86400;%conversion to mm/day
312			cc=[[0:5:25] 35 [50:25:200] [250 300]]*0.04; title0='std(EMPMR to ocean+ice)';
313			if doAnomalies; cc=scaleAnom[0:0.1:1]0.5; end;
314			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
315			myCaption={myYmeanTxt,' standard deviation -- EMPMR to ocean+ice (W/m2)'};
316			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
317
318			%empmr from ocean:
319			fld=alldiag.oceFWflx_std*86400;%conversion to mm/day
320			cc=[[0:5:25] 35 [50:25:200] [250 300]]*0.04; title0='std(EMPMR to ocean)';
321			if doAnomalies; cc=scaleAnom[0:0.1:1]0.5; end;
322			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
323			myCaption={myYmeanTxt,' standard deviation -- EMPMR to ocean (W/m2)'};
324			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
325			end;
326
327			end;
328
329			if sum(strcmp(choicePlot,'all'))\|sum(strcmp(choicePlot,'tau'));
330
331	gforget	1.6	if addToTex; write2tex(fileTex,1,'surface wind stress',2); end;
332
333	gforget	1.4	%zonal wind stress:
334			fld=alldiag.fldTZ_mean;
335			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='zonal wind stress';
336			if doAnomalies; cc=scaleAnom[-1:0.1:1]0.01; end;
337			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
338			myCaption={myYmeanTxt,'mean -- zonal wind stress (N/m2)'};
339			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
340
341			%meridional wind stress:
342			fld=alldiag.fldTM_mean;
343			cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='meridional wind stress';
344			if doAnomalies; cc=scaleAnom[-1:0.1:1]0.01; end;
345			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
346			myCaption={myYmeanTxt,'mean -- meridional wind stress (N/m2)'};
347			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
348
349	gforget	1.6	%fld=alldiag.curlTau_mean;
350			%cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/5e8; title0='wind stress curl';
351			%if doAnomalies; cc=scaleAnom[-1:0.1:1]0.01; end;
352			%figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
353			%myCaption={myYmeanTxt,'mean -- wind stress curl (N/m3)'};
354			%if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
355	gforget	1.4
356			if multiTimes;
357			%zonal wind stress:
358			fld=alldiag.fldTZ_std;
359			cc=[[0:5:25] 35 [50:25:200] [250 300]]/2000; title0='std(tauZ)';
360			if doAnomalies; cc=scaleAnom[0:0.1:1]0.005; end;
361			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
362			myCaption={myYmeanTxt,' standard deviation -- tauZ (W/m2)'};
363			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
364
365			%meridional wind stress:
366			fld=alldiag.fldTM_std;
367			cc=[[0:5:25] 35 [50:25:200] [250 300]]/2000; title0='std(tauM)';
368			if doAnomalies; cc=scaleAnom[0:0.1:1]0.005; end;
369			figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
370			myCaption={myYmeanTxt,' standard deviation -- tauM (W/m2)'};
371			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
372
373	gforget	1.6	%fld=alldiag.curlTau_std;
374			%cc=[[0:5:25] 35 [50:25:200] [250 300]]/1e9; title0='wind stress curl';
375			%if doAnomalies; cc=scaleAnom[-1:0.1:1]0.01; end;
376			%figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
377			%myCaption={myYmeanTxt,'standard deviation -- tauCurl (N/m3)'};
378			%if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
379	gforget	1.4	end;
380
381			end;
382
383	gforget	1.1	end;
384	gforget	1.4
385