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 SIatmQnt_2s SIatmFW_s SIatmFW_2s oceQnet_s oceQnet_2s ' ...
               'oceFWflx_s oceFWflx_2s fldTZ_s fldTZ_2s fldTM_s fldTM_2s ' ...
               'curlTau_s curlTau_2s fldETAN_s fldETAN_2s ' ...
               'fldSSH_s fldSSH_2s fldMLD_s fldMLD_2s'];
    if ~doOmit3Dfields;
      listDiags=[listDiags ' THETA_s THETA_2s SALT_s SALT_2s WVELMASS_s 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/' ],[dirModel 'diags/']};

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);
    fldSSH=(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>myparms.recInAve(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{1+2*(jj-1)}(1:end-2);
        if ii==myparms.recInAve(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

  nfiles=length(dir(fullfile(dirMat,[fileMat '_*.mat'])));
  if nfiles>=myparms.recInAve(2);  
    ifile1=myparms.recInAve(2);
    ifile0=myparms.recInAve(1)-1;
  else;
    ifile1=nfiles;
    ifile0=0;
  end;

  %load last cumsum
  alldiag=diags_read_from_mat(dirMat,[fileMat '_*.mat'],'',ifile1);
  %load first cumsum if needed
  if ifile0>0;
    tmpdiag=diags_read_from_mat(dirMat,[fileMat '_*.mat'],'',ifile0);
    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;
  %ensure backward compatibility
  if ~doOmit3Dfields;
    if ~isfield(alldiag,'WVELMASS_s'); alldiag.WVELMASS_s=NaN*alldiag.THETA_s; end;
    if ~isfield(alldiag,'WVELMASS_2s'); alldiag.WVELMASS_2s=NaN*alldiag.THETA_2s; end;
  end;
  if ~isfield(alldiag,'fldSSH_s'); alldiag.fldSSH_s=alldiag.fldETANLEADS_s; end;
  if ~isfield(alldiag,'fldSSH_2s'); alldiag.fldSSH_2s=alldiag.fldETANLEADS_2s; end;
  jj=find(strcmp(alldiag.listDiags,'fldETANLEADS_s'));
  if ~isempty(jj); alldiag.listDiags{jj}='fldSSH_s'; end;
  jj=find(strcmp(alldiag.listDiags,'fldETANLEADS_2s'));
  if ~isempty(jj); alldiag.listDiags{jj}='fldSSH_2s'; end;
  %
  n=diff(myparms.recInAve)+1;
  for jj=1:length(alldiag.listDiags)/2;
     tmp0=alldiag.listDiags{1+2*(jj-1)};
     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.fldSSH_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.fldSSH_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 degC)'; 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 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/m$^2$)'};
    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/m$^2$)'};
    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/m$^2$)'};
        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/m$^2$)'};
        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 freshwater 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='E-P-R 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 -- E-P-R 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='E-P-R 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 -- E-P-R 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(E-P-R 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 -- E-P-R to ocean+ice (W/m$^2$)'};
        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(E-P-R 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 -- E-P-R to ocean (W/m$^2$)'};
        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/m$^2$)'};
    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/m$^2$)'};
    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/m$^3$)'};
    %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/m$^2$)'};
        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/m$^2$)'};
        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/m$^3$)'};
        %if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    end;

    end;

end;


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