matlab_class/gcmfaces_diags/diags_set_LAYERS.m


layersOffline=1;
layersEulerian=0;
integrateFromBottom=1;
if isempty(whos('setDiagsParams')); setDiagsParams={}; end;

if length(setDiagsParams)==2;
    layersName=setDiagsParams{1};
    layersLims=setDiagsParams{2};
    layersOffline=1;
    layersEulerian=1;%hacked : should be 0
elseif length(setDiagsParams)==1;
    layersName=setDiagsParams{1};
else;
    layersName='sigma';
    pref0=0;
    layersLims0=[18:0.5:22.5 22.8:0.2:27.4 27.5:0.02:28 28.2:0.2:29.2 29.5:0.5:31]';
    %pref0=2000;
    %layersLims0=[30:0.1:38]';
end;

layersFileSuff=['LAYERS_' layersName];
if layersOffline; layersFileSuff=['LAYERS_offline_' layersName]; end;

%override default file name:
%---------------------------
fileMat=['diags_set_' layersFileSuff];

if userStep==1;%diags to be computed
    listDiags=['layersParams gloOV gloThick gloDelThick gloBAR gloMT gloD'];
    listBasins=1;
    if sum([90 1170]~=mygrid.ioSize)==0;
        listBasins=[1:3];
        mskC=v4_basin({'atlExt'});
        if isempty(mskC); listBasins=1; end;
    end;
    if length(listBasins)==3;
      listDiags=[listDiags ' atlOV atlThick atlDelThick atlBAR atlMT atlD'];
      listDiags=[listDiags ' pacindOV pacindThick pacindDelThick pacindBAR pacindMT pacindD'];
    end;
elseif userStep==2&~layersOffline;%input files and variables
    listFlds={    ['LaUH' layersName],['LaVH' layersName],...
        ['LaHw' layersName],['LaHs' layersName]};
    listFldsNames=deblank(listFlds);
    listFiles={'layersDiags'};
    listSubdirs={[dirModel 'diags/LAYERS/' ],[dirModel 'diags/']};
    %load layersLims consistent with online MITgcmpkg/layers
    layersLims=squeeze(rdmds([dirModel 'diags/LAYERS/layers' layersName]));
elseif userStep==2&layersOffline;%input files and variables
    listFlds={'THETA','SALT','UVELMASS','VVELMASS','GM_PsiX','GM_PsiY'};
    listFldsNames=deblank(listFlds);
    listFiles={'state_3d_set1','trsp_3d_set1','trsp_3d_set2'};
elseif userStep==3;

    %override default file name:
    %---------------------------
    fileMat=['diags_set_' layersFileSuff '_' num2str(tt) '.mat'];
    
    if ~layersOffline;
        
        eval(['fldU=LaUH' layersName '; fldV=LaVH' layersName ';']);
        if ~isempty(whos(['LaHw' layersName]));
            eval(['fldHw=LaHw' layersName '; fldHs=LaHs' layersName ';']);
        else;
            fldHw=[]; fldHs=[];
        end;
        fldD=[];
        layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
        
    else;
        
        if ~layersEulerian;
            [fldUbolus,fldVbolus,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY);
            fldUbolus=fldUbolus.*mygrid.mskW; fldVbolus=fldVbolus.*mygrid.mskS;
            UVELMASS=UVELMASS+fldUbolus; VVELMASS=VVELMASS+fldVbolus;
        end;
        
        U=UVELMASS.*mk3D(mygrid.DRF,UVELMASS);
        V=VVELMASS.*mk3D(mygrid.DRF,VVELMASS);

        if strcmp(layersName,'theta');
            tracer=THETA;
            if isempty(whos('layersLims')); layersLims=[-2:35]; end;
        elseif strcmp(layersName,'salt');
            tracer=SALT;
            if isempty(whos('layersLims')); layersLims=[33:0.1:38]; end;
        elseif strcmp(layersName,'sigma');
            P=mk3D(-mygrid.RC,THETA);
            t=convert2vector(THETA);
            s=convert2vector(SALT);
            p=convert2vector(P);
            pref=pref0+0*p;
            [rhop,rhpis,rhor] = density(t(:),s(:),p(:),pref(:));
            rhor=reshape(rhor,size(t));
            tracer=convert2vector(rhor)-1000;
            if isempty(whos('layersLims')); layersLims=layersLims0; end;
        end
        
        layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
        [fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,2);
        fldHw=[]; fldHs=[];
        [fldD]=layers_remap(P,'intensive',tracer,layersGrid,2);

    end;%if ~layersOffline;

    for bb=listBasins;

    %mask : global, atlantic or Pac+Ind
    if bb==1;       mskC=mygrid.mskC(:,:,1); mskW=mygrid.mskW(:,:,1); mskS=mygrid.mskS(:,:,1);
    elseif bb==2;   [mskC,mskW,mskS]=v4_basin({'atlExt'});
    elseif bb==3;   [mskC,mskW,mskS]=v4_basin({'pacExt','indExt'});
    end;
    mskC3d=1*(mskC>0); mskW3d=1*(mskW>0); mskS3d=1*(mskS>0);
    mskC3d=mk3D(mskC3d,fldU); mskW3d=mk3D(mskW3d,fldU); mskS3d=mk3D(mskS3d,fldU);

    %the overturning streamfunction computation itself
    layersOV=calc_overturn(fldU.*mskW3d,fldV.*mskS3d,1,{'dh'});

    %the associated barotropic streamfunction (for checking)
    layersBAR=calc_barostream(fldU.*mskW3d,fldV.*mskS3d,1,{'dh'});

    %meridional transport per layer:
    layersMT=diff(layersOV,1,2); 

    if ~isempty(fldD);
      %compute zonal mean depth from fldD
      layersD=calc_zonmean_T(fldD.*mskC3d);
    else;
      layersD=[];
    end;

    %the associated thickness
    if isempty(fldHw);
        layersThick=[];
        layersDelThick=[];
    else;
        %interpolate to tracer points
        fldH=0*fldHw;
        [fldHwp,fldHsp]=exch_UV(LaHw1SLT.*mskW3d,LaHs1SLT.*mskS3d);
        for iF=1:fldH.nFaces;
            tmpA=fldHwp{iF}(2:end,:,:);
            tmpB=fldHwp{iF}(1:end-1,:,:);
            tmpC=fldHsp{iF}(:,2:end,:);
            tmpD=fldHsp{iF}(:,1:end-1,:);
            tmpTot=tmpA+tmpB+tmpC+tmpD;
            tmpNb=1*(tmpA~=0)+1*(tmpB~=0)+1*(tmpC~=0)+1*(tmpD~=0);
            jj=find(tmpNb>0); tmpTot(jj)=tmpTot(jj)./tmpNb(jj);
            jj=find(isnan(tmpTot)); tmpTot(jj)=0;
            fldH{iF}=tmpTot;
        end;
        %compute zonal mean
        fldH=calc_zonmean_T(fldH.*mskC3d);
        %integrate to overturning points
        if integrateFromBottom;
            layersThick=[flipdim(cumsum(flipdim(fldH,2),2),2) zeros(size(fldH,1),1)];
        else;
            layersThick=[zeros(size(fldH,1),1) cumsum(fldH,2)];
        end;
        %compute dH/dLayer
        layersDelThick=[zeros(size(fldH,1),1) fldH./( ones(size(fldH,1),1)*diff(layersLims') )];
    end;

    %store to global, atlantic or Pac+Ind arrays:
    if bb==1;
        gloOV=layersOV; gloThick=layersThick; gloDelThick=layersDelThick; 
        gloBAR=layersBAR; gloMT=layersMT; gloD=layersD;
    elseif bb==2;
        kk=find(mygrid.LATS<-35|mygrid.LATS>70);
        layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
        if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
        atlOV=layersOV; atlThick=layersThick; atlDelThick=layersDelThick;
        atlBAR=layersBAR; atlMT=layersMT; atlD=layersD;
    elseif bb==3;
        kk=find(mygrid.LATS<-35|mygrid.LATS>65);
        layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
        if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
        pacindOV=layersOV; pacindThick=layersThick; pacindDelThick=layersDelThick;
        pacindBAR=layersBAR; pacindMT=layersMT; pacindD=layersD;
    end;
    end;
    
    layersParams.name=layersName;
    layersParams.LC=layersGrid;
    layersParams.LF=layersLims;
    layersParams.isOffline=layersOffline;
    layersParams.suffix=layersFileSuff;
    layersParams.isEulerian=layersEulerian;

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

elseif userStep==-1;%plotting
   
    if ~sum(strcmp(listDiags,'atlOV')); multiBasins=0; end;
 
    X=mygrid.LATS*ones(1,length(alldiag.layersParams(1).LF));
    Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams(1).LF');
    cc=[[-50:10:-30] [-24:3:24] [30:10:50]];

    %meridional streamfunction (Eulerian) :
    fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN; title0='Meridional Stream Function';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
    figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
    pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
    myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    if multiBasins;

    %meridional streamfunction in Atlantic:
    fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN; title0='Atlantic Meridional Stream Function';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
    figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
    pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
    myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    %meridional streamfunction Pacific+Indian:
    fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN; title0='Pac+Ind Meridional Stream Function';
    if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
    figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
    pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
    myCaption={myYmeanTxt,'mean -- Pac+Ind overturning streamfunction (Sv)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    end;

end;
1	gforget	1.1
2	gforget	1.5	layersOffline=1;
3	gforget	1.4	layersEulerian=0;
4	gforget	1.1	integrateFromBottom=1;
5			if isempty(whos('setDiagsParams')); setDiagsParams={}; end;
6
7	gforget	1.3	if length(setDiagsParams)==2;
8	gforget	1.1	layersName=setDiagsParams{1};
9	gforget	1.3	layersLims=setDiagsParams{2};
10	gforget	1.1	layersOffline=1;
11	gforget	1.4	layersEulerian=1;%hacked : should be 0
12	gforget	1.1	elseif length(setDiagsParams)==1;
13	gforget	1.3	layersName=setDiagsParams{1};
14	gforget	1.1	else;
15	gforget	1.8	layersName='sigma';
16	gforget	1.9	pref0=0;
17			layersLims0=[18:0.5:22.5 22.8:0.2:27.4 27.5:0.02:28 28.2:0.2:29.2 29.5:0.5:31]';
18			%pref0=2000;
19			%layersLims0=[30:0.1:38]';
20	gforget	1.1	end;
21
22			layersFileSuff=['LAYERS_' layersName];
23			if layersOffline; layersFileSuff=['LAYERS_offline_' layersName]; end;
24	gforget	1.3
25			%override default file name:
26			%---------------------------
27			fileMat=['diags_set_' layersFileSuff];
28	gforget	1.1
29			if userStep==1;%diags to be computed
30	gforget	1.5	listDiags=['layersParams gloOV gloThick gloDelThick gloBAR gloMT gloD'];
31	gforget	1.8	listBasins=1;
32	gforget	1.3	if sum([90 1170]~=mygrid.ioSize)==0;
33	gforget	1.8	listBasins=[1:3];
34			mskC=v4_basin({'atlExt'});
35			if isempty(mskC); listBasins=1; end;
36			end;
37			if length(listBasins)==3;
38	gforget	1.5	listDiags=[listDiags ' atlOV atlThick atlDelThick atlBAR atlMT atlD'];
39			listDiags=[listDiags ' pacindOV pacindThick pacindDelThick pacindBAR pacindMT pacindD'];
40	gforget	1.3	end;
41	gforget	1.1	elseif userStep==2&~layersOffline;%input files and variables
42			listFlds={ ['LaUH' layersName],['LaVH' layersName],...
43			['LaHw' layersName],['LaHs' layersName]};
44			listFldsNames=deblank(listFlds);
45			listFiles={'layersDiags'};
46	gforget	1.7	listSubdirs={[dirModel 'diags/LAYERS/' ],[dirModel 'diags/']};
47	gforget	1.1	%load layersLims consistent with online MITgcmpkg/layers
48			layersLims=squeeze(rdmds([dirModel 'diags/LAYERS/layers' layersName]));
49			elseif userStep==2&layersOffline;%input files and variables
50			listFlds={'THETA','SALT','UVELMASS','VVELMASS','GM_PsiX','GM_PsiY'};
51			listFldsNames=deblank(listFlds);
52			listFiles={'state_3d_set1','trsp_3d_set1','trsp_3d_set2'};
53	gforget	1.3	elseif userStep==3;
54	gforget	1.1
55			%override default file name:
56			%---------------------------
57	gforget	1.3	fileMat=['diags_set_' layersFileSuff '_' num2str(tt) '.mat'];
58	gforget	1.1
59			if ~layersOffline;
60
61			eval(['fldU=LaUH' layersName '; fldV=LaVH' layersName ';']);
62			if ~isempty(whos(['LaHw' layersName]));
63			eval(['fldHw=LaHw' layersName '; fldHs=LaHs' layersName ';']);
64			else;
65			fldHw=[]; fldHs=[];
66			end;
67	gforget	1.5	fldD=[];
68	gforget	1.1	layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
69
70			else;
71
72			if ~layersEulerian;
73			[fldUbolus,fldVbolus,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY);
74			fldUbolus=fldUbolus.mygrid.mskW; fldVbolus=fldVbolus.mygrid.mskS;
75			UVELMASS=UVELMASS+fldUbolus; VVELMASS=VVELMASS+fldVbolus;
76			end;
77
78			U=UVELMASS.*mk3D(mygrid.DRF,UVELMASS);
79			V=VVELMASS.*mk3D(mygrid.DRF,VVELMASS);
80	gforget	1.8
81	gforget	1.5	if strcmp(layersName,'theta');
82	gforget	1.1	tracer=THETA;
83	gforget	1.4	if isempty(whos('layersLims')); layersLims=[-2:35]; end;
84	gforget	1.5	elseif strcmp(layersName,'salt');
85	gforget	1.1	tracer=SALT;
86	gforget	1.4	if isempty(whos('layersLims')); layersLims=[33:0.1:38]; end;
87	gforget	1.8	elseif strcmp(layersName,'sigma');
88	gforget	1.1	P=mk3D(-mygrid.RC,THETA);
89			t=convert2vector(THETA);
90			s=convert2vector(SALT);
91			p=convert2vector(P);
92	gforget	1.8	pref=pref0+0*p;
93	gforget	1.10	[rhop,rhpis,rhor] = density(t(:),s(:),p(:),pref(:));
94			rhor=reshape(rhor,size(t));
95	gforget	1.6	tracer=convert2vector(rhor)-1000;
96	gforget	1.8	if isempty(whos('layersLims')); layersLims=layersLims0; end;
97	gforget	1.1	end
98
99			layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
100	gforget	1.5	[fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,2);
101	gforget	1.1	fldHw=[]; fldHs=[];
102	gforget	1.6	[fldD]=layers_remap(P,'intensive',tracer,layersGrid,2);
103	gforget	1.5
104	gforget	1.1	end;%if ~layersOffline;
105	gforget	1.3
106			for bb=listBasins;
107
108			%mask : global, atlantic or Pac+Ind
109			if bb==1; mskC=mygrid.mskC(:,:,1); mskW=mygrid.mskW(:,:,1); mskS=mygrid.mskS(:,:,1);
110			elseif bb==2; [mskC,mskW,mskS]=v4_basin({'atlExt'});
111			elseif bb==3; [mskC,mskW,mskS]=v4_basin({'pacExt','indExt'});
112	gforget	1.1	end;
113			mskC3d=1(mskC>0); mskW3d=1(mskW>0); mskS3d=1*(mskS>0);
114			mskC3d=mk3D(mskC3d,fldU); mskW3d=mk3D(mskW3d,fldU); mskS3d=mk3D(mskS3d,fldU);
115	gforget	1.2
116	gforget	1.3	%the overturning streamfunction computation itself
117			layersOV=calc_overturn(fldU.mskW3d,fldV.mskS3d,1,{'dh'});
118
119	gforget	1.1	%the associated barotropic streamfunction (for checking)
120	gforget	1.3	layersBAR=calc_barostream(fldU.mskW3d,fldV.mskS3d,1,{'dh'});
121	gforget	1.5
122			%meridional transport per layer:
123	gforget	1.6	layersMT=diff(layersOV,1,2);
124	gforget	1.5
125	gforget	1.6	if ~isempty(fldD);
126			%compute zonal mean depth from fldD
127			layersD=calc_zonmean_T(fldD.*mskC3d);
128	gforget	1.5	else;
129			layersD=[];
130			end;
131
132	gforget	1.1	%the associated thickness
133			if isempty(fldHw);
134			layersThick=[];
135			layersDelThick=[];
136			else;
137			%interpolate to tracer points
138			fldH=0*fldHw;
139	gforget	1.3	[fldHwp,fldHsp]=exch_UV(LaHw1SLT.mskW3d,LaHs1SLT.mskS3d);
140	gforget	1.1	for iF=1:fldH.nFaces;
141			tmpA=fldHwp{iF}(2:end,:,:);
142			tmpB=fldHwp{iF}(1:end-1,:,:);
143			tmpC=fldHsp{iF}(:,2:end,:);
144			tmpD=fldHsp{iF}(:,1:end-1,:);
145			tmpTot=tmpA+tmpB+tmpC+tmpD;
146			tmpNb=1(tmpA~=0)+1(tmpB~=0)+1(tmpC~=0)+1(tmpD~=0);
147			jj=find(tmpNb>0); tmpTot(jj)=tmpTot(jj)./tmpNb(jj);
148			jj=find(isnan(tmpTot)); tmpTot(jj)=0;
149			fldH{iF}=tmpTot;
150			end;
151			%compute zonal mean
152			fldH=calc_zonmean_T(fldH.*mskC3d);
153			%integrate to overturning points
154			if integrateFromBottom;
155			layersThick=[flipdim(cumsum(flipdim(fldH,2),2),2) zeros(size(fldH,1),1)];
156			else;
157			layersThick=[zeros(size(fldH,1),1) cumsum(fldH,2)];
158			end;
159			%compute dH/dLayer
160			layersDelThick=[zeros(size(fldH,1),1) fldH./( ones(size(fldH,1),1)*diff(layersLims') )];
161			end;
162	gforget	1.3
163			%store to global, atlantic or Pac+Ind arrays:
164			if bb==1;
165	gforget	1.5	gloOV=layersOV; gloThick=layersThick; gloDelThick=layersDelThick;
166	gforget	1.6	gloBAR=layersBAR; gloMT=layersMT; gloD=layersD;
167	gforget	1.3	elseif bb==2;
168			kk=find(mygrid.LATS<-35\|mygrid.LATS>70);
169	gforget	1.5	layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
170			if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
171			atlOV=layersOV; atlThick=layersThick; atlDelThick=layersDelThick;
172			atlBAR=layersBAR; atlMT=layersMT; atlD=layersD;
173	gforget	1.3	elseif bb==3;
174			kk=find(mygrid.LATS<-35\|mygrid.LATS>65);
175	gforget	1.5	layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
176			if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
177			pacindOV=layersOV; pacindThick=layersThick; pacindDelThick=layersDelThick;
178			pacindBAR=layersBAR; pacindMT=layersMT; pacindD=layersD;
179	gforget	1.3	end;
180			end;
181	gforget	1.1
182			layersParams.name=layersName;
183			layersParams.LC=layersGrid;
184			layersParams.LF=layersLims;
185			layersParams.isOffline=layersOffline;
186			layersParams.suffix=layersFileSuff;
187			layersParams.isEulerian=layersEulerian;
188	gforget	1.3
189			%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
190			%===================== PLOTTING SEQUENCE BEGINS =========================%
191
192			elseif userStep==-1;%plotting
193	gforget	1.8
194			if ~sum(strcmp(listDiags,'atlOV')); multiBasins=0; end;
195
196			X=mygrid.LATS*ones(1,length(alldiag.layersParams(1).LF));
197			Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams(1).LF');
198			cc=[[-50:10:-30] [-24:3:24] [30:10:50]];
199	gforget	1.3
200			%meridional streamfunction (Eulerian) :
201	gforget	1.8	fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN; title0='Meridional Stream Function';
202	gforget	1.3	if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
203			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
204	gforget	1.8	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
205	gforget	1.3	myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'};
206			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
207
208			if multiBasins;
209	gforget	1.8
210			%meridional streamfunction in Atlantic:
211			fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN; title0='Atlantic Meridional Stream Function';
212	gforget	1.3	if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
213			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
214	gforget	1.8	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
215	gforget	1.3	myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'};
216			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
217	gforget	1.8
218			%meridional streamfunction Pacific+Indian:
219			fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN; title0='Pac+Ind Meridional Stream Function';
220	gforget	1.3	if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
221			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
222	gforget	1.8	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
223	gforget	1.3	myCaption={myYmeanTxt,'mean -- Pac+Ind overturning streamfunction (Sv)'};
224			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
225	gforget	1.8
226	gforget	1.3	end;
227
228	gforget	1.1	end;