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='1SLT';
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'];
    if sum([90 1170]~=mygrid.ioSize)==0;
      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/' ]};
    %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,'sig0');
            P=mk3D(-mygrid.RC,THETA);
            t=convert2vector(THETA);
            s=convert2vector(SALT);
            p=convert2vector(P);
            %pref=2000+0*p;%could be made an extra param
            pref=0*p;%could be made an extra param
            [rhop,rhpis,rhor] = density(t,s,p,pref);
            tracer=convert2vector(rhor)-1000;
            %if isempty(whos('layersLims')); layersLims=1000+[30:0.1:38]'; end;
            if isempty(whos('layersLims')); 
              layersLims=[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]'; 
            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;

    if sum([90 1170]~=mygrid.ioSize)==0;
        listBasins=[1:3];
    else;
        listBasins=1;
    end;    

    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
    

    %meridional streamfunction (Eulerian) :
    fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN;
    X=mygrid.LATS*ones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams.LF');
    cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; 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);
    myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

    if multiBasins;
    %meridional streamfunction (Eulerian):
    fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN;
    X=mygrid.LATS*ones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams.LF');
    cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; 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);
    myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'};
    if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
    %meridional streamfunction (residual):
    fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN;
    X=mygrid.LATS*ones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams.LF');
    cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; 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);
    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.3	layersName='1SLT';
16	gforget	1.1	end;
17
18			layersFileSuff=['LAYERS_' layersName];
19			if layersOffline; layersFileSuff=['LAYERS_offline_' layersName]; end;
20	gforget	1.3
21			%override default file name:
22			%---------------------------
23			fileMat=['diags_set_' layersFileSuff];
24	gforget	1.1
25			if userStep==1;%diags to be computed
26	gforget	1.5	listDiags=['layersParams gloOV gloThick gloDelThick gloBAR gloMT gloD'];
27	gforget	1.3	if sum([90 1170]~=mygrid.ioSize)==0;
28	gforget	1.5	listDiags=[listDiags ' atlOV atlThick atlDelThick atlBAR atlMT atlD'];
29			listDiags=[listDiags ' pacindOV pacindThick pacindDelThick pacindBAR pacindMT pacindD'];
30	gforget	1.3	end;
31	gforget	1.1	elseif userStep==2&~layersOffline;%input files and variables
32			listFlds={ ['LaUH' layersName],['LaVH' layersName],...
33			['LaHw' layersName],['LaHs' layersName]};
34			listFldsNames=deblank(listFlds);
35			listFiles={'layersDiags'};
36			listSubdirs={[dirModel 'diags/LAYERS/' ]};
37			%load layersLims consistent with online MITgcmpkg/layers
38			layersLims=squeeze(rdmds([dirModel 'diags/LAYERS/layers' layersName]));
39			elseif userStep==2&layersOffline;%input files and variables
40			listFlds={'THETA','SALT','UVELMASS','VVELMASS','GM_PsiX','GM_PsiY'};
41			listFldsNames=deblank(listFlds);
42			listFiles={'state_3d_set1','trsp_3d_set1','trsp_3d_set2'};
43	gforget	1.3	elseif userStep==3;
44	gforget	1.1
45			%override default file name:
46			%---------------------------
47	gforget	1.3	fileMat=['diags_set_' layersFileSuff '_' num2str(tt) '.mat'];
48	gforget	1.1
49			if ~layersOffline;
50
51			eval(['fldU=LaUH' layersName '; fldV=LaVH' layersName ';']);
52			if ~isempty(whos(['LaHw' layersName]));
53			eval(['fldHw=LaHw' layersName '; fldHs=LaHs' layersName ';']);
54			else;
55			fldHw=[]; fldHs=[];
56			end;
57	gforget	1.5	fldD=[];
58	gforget	1.1	layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
59
60			else;
61
62			if ~layersEulerian;
63			[fldUbolus,fldVbolus,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY);
64			fldUbolus=fldUbolus.mygrid.mskW; fldVbolus=fldVbolus.mygrid.mskS;
65			UVELMASS=UVELMASS+fldUbolus; VVELMASS=VVELMASS+fldVbolus;
66			end;
67
68			U=UVELMASS.*mk3D(mygrid.DRF,UVELMASS);
69			V=VVELMASS.*mk3D(mygrid.DRF,VVELMASS);
70
71	gforget	1.5	if strcmp(layersName,'theta');
72	gforget	1.1	tracer=THETA;
73	gforget	1.4	if isempty(whos('layersLims')); layersLims=[-2:35]; end;
74	gforget	1.5	elseif strcmp(layersName,'salt');
75	gforget	1.1	tracer=SALT;
76	gforget	1.4	if isempty(whos('layersLims')); layersLims=[33:0.1:38]; end;
77	gforget	1.5	elseif strcmp(layersName,'sig0');
78	gforget	1.1	P=mk3D(-mygrid.RC,THETA);
79			t=convert2vector(THETA);
80			s=convert2vector(SALT);
81			p=convert2vector(P);
82	gforget	1.5	%pref=2000+0*p;%could be made an extra param
83			pref=0*p;%could be made an extra param
84	gforget	1.1	[rhop,rhpis,rhor] = density(t,s,p,pref);
85	gforget	1.6	tracer=convert2vector(rhor)-1000;
86	gforget	1.5	%if isempty(whos('layersLims')); layersLims=1000+[30:0.1:38]'; end;
87	gforget	1.6	if isempty(whos('layersLims'));
88			layersLims=[18:0.5:22.5 22.8:0.2:27.4 ...
89			27.5:0.02:28 28.2:0.2:29.2 29.5:0.5:31]';
90			end;
91	gforget	1.1	end
92
93			layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
94	gforget	1.5	[fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,2);
95	gforget	1.1	fldHw=[]; fldHs=[];
96	gforget	1.6	[fldD]=layers_remap(P,'intensive',tracer,layersGrid,2);
97	gforget	1.5
98	gforget	1.1	end;%if ~layersOffline;
99	gforget	1.3
100			if sum([90 1170]~=mygrid.ioSize)==0;
101			listBasins=[1:3];
102			else;
103			listBasins=1;
104			end;
105
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
194
195			%meridional streamfunction (Eulerian) :
196			fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN;
197			X=mygrid.LATSones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)(alldiag.layersParams.LF');
198			cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Meridional Stream Function';
199			if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
200			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
201			pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0);
202			myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'};
203			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
204
205			if multiBasins;
206			%meridional streamfunction (Eulerian):
207			fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN;
208			X=mygrid.LATSones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)(alldiag.layersParams.LF');
209			cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Atlantic Meridional Stream Function';
210			if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
211			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
212			pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0);
213			myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'};
214			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
215			%meridional streamfunction (residual):
216			fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN;
217			X=mygrid.LATSones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)(alldiag.layersParams.LF');
218			cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Pac+Ind Meridional Stream Function';
219			if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
220			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
221			pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0);
222			myCaption={myYmeanTxt,'mean -- Pac+Ind overturning streamfunction (Sv)'};
223			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
224			end;
225
226	gforget	1.1
227			end;