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);
            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			%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.1	[rhop,rhpis,rhor] = density(t,s,p,pref);
94	gforget	1.6	tracer=convert2vector(rhor)-1000;
95	gforget	1.8	if isempty(whos('layersLims')); layersLims=layersLims0; end;
96	gforget	1.1	end
97
98			layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
99	gforget	1.5	[fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,2);
100	gforget	1.1	fldHw=[]; fldHs=[];
101	gforget	1.6	[fldD]=layers_remap(P,'intensive',tracer,layersGrid,2);
102	gforget	1.5
103	gforget	1.1	end;%if ~layersOffline;
104	gforget	1.3
105			for bb=listBasins;
106
107			%mask : global, atlantic or Pac+Ind
108			if bb==1; mskC=mygrid.mskC(:,:,1); mskW=mygrid.mskW(:,:,1); mskS=mygrid.mskS(:,:,1);
109			elseif bb==2; [mskC,mskW,mskS]=v4_basin({'atlExt'});
110			elseif bb==3; [mskC,mskW,mskS]=v4_basin({'pacExt','indExt'});
111	gforget	1.1	end;
112			mskC3d=1(mskC>0); mskW3d=1(mskW>0); mskS3d=1*(mskS>0);
113			mskC3d=mk3D(mskC3d,fldU); mskW3d=mk3D(mskW3d,fldU); mskS3d=mk3D(mskS3d,fldU);
114	gforget	1.2
115	gforget	1.3	%the overturning streamfunction computation itself
116			layersOV=calc_overturn(fldU.mskW3d,fldV.mskS3d,1,{'dh'});
117
118	gforget	1.1	%the associated barotropic streamfunction (for checking)
119	gforget	1.3	layersBAR=calc_barostream(fldU.mskW3d,fldV.mskS3d,1,{'dh'});
120	gforget	1.5
121			%meridional transport per layer:
122	gforget	1.6	layersMT=diff(layersOV,1,2);
123	gforget	1.5
124	gforget	1.6	if ~isempty(fldD);
125			%compute zonal mean depth from fldD
126			layersD=calc_zonmean_T(fldD.*mskC3d);
127	gforget	1.5	else;
128			layersD=[];
129			end;
130
131	gforget	1.1	%the associated thickness
132			if isempty(fldHw);
133			layersThick=[];
134			layersDelThick=[];
135			else;
136			%interpolate to tracer points
137			fldH=0*fldHw;
138	gforget	1.3	[fldHwp,fldHsp]=exch_UV(LaHw1SLT.mskW3d,LaHs1SLT.mskS3d);
139	gforget	1.1	for iF=1:fldH.nFaces;
140			tmpA=fldHwp{iF}(2:end,:,:);
141			tmpB=fldHwp{iF}(1:end-1,:,:);
142			tmpC=fldHsp{iF}(:,2:end,:);
143			tmpD=fldHsp{iF}(:,1:end-1,:);
144			tmpTot=tmpA+tmpB+tmpC+tmpD;
145			tmpNb=1(tmpA~=0)+1(tmpB~=0)+1(tmpC~=0)+1(tmpD~=0);
146			jj=find(tmpNb>0); tmpTot(jj)=tmpTot(jj)./tmpNb(jj);
147			jj=find(isnan(tmpTot)); tmpTot(jj)=0;
148			fldH{iF}=tmpTot;
149			end;
150			%compute zonal mean
151			fldH=calc_zonmean_T(fldH.*mskC3d);
152			%integrate to overturning points
153			if integrateFromBottom;
154			layersThick=[flipdim(cumsum(flipdim(fldH,2),2),2) zeros(size(fldH,1),1)];
155			else;
156			layersThick=[zeros(size(fldH,1),1) cumsum(fldH,2)];
157			end;
158			%compute dH/dLayer
159			layersDelThick=[zeros(size(fldH,1),1) fldH./( ones(size(fldH,1),1)*diff(layersLims') )];
160			end;
161	gforget	1.3
162			%store to global, atlantic or Pac+Ind arrays:
163			if bb==1;
164	gforget	1.5	gloOV=layersOV; gloThick=layersThick; gloDelThick=layersDelThick;
165	gforget	1.6	gloBAR=layersBAR; gloMT=layersMT; gloD=layersD;
166	gforget	1.3	elseif bb==2;
167			kk=find(mygrid.LATS<-35\|mygrid.LATS>70);
168	gforget	1.5	layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
169			if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
170			atlOV=layersOV; atlThick=layersThick; atlDelThick=layersDelThick;
171			atlBAR=layersBAR; atlMT=layersMT; atlD=layersD;
172	gforget	1.3	elseif bb==3;
173			kk=find(mygrid.LATS<-35\|mygrid.LATS>65);
174	gforget	1.5	layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
175			if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
176			pacindOV=layersOV; pacindThick=layersThick; pacindDelThick=layersDelThick;
177			pacindBAR=layersBAR; pacindMT=layersMT; pacindD=layersD;
178	gforget	1.3	end;
179			end;
180	gforget	1.1
181			layersParams.name=layersName;
182			layersParams.LC=layersGrid;
183			layersParams.LF=layersLims;
184			layersParams.isOffline=layersOffline;
185			layersParams.suffix=layersFileSuff;
186			layersParams.isEulerian=layersEulerian;
187	gforget	1.3
188			%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
189			%===================== PLOTTING SEQUENCE BEGINS =========================%
190
191			elseif userStep==-1;%plotting
192	gforget	1.8
193			if ~sum(strcmp(listDiags,'atlOV')); multiBasins=0; end;
194
195			X=mygrid.LATS*ones(1,length(alldiag.layersParams(1).LF));
196			Y=ones(length(mygrid.LATS),1)*(alldiag.layersParams(1).LF');
197			cc=[[-50:10:-30] [-24:3:24] [30:10:50]];
198	gforget	1.3
199			%meridional streamfunction (Eulerian) :
200	gforget	1.8	fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN; title0='Meridional Stream Function';
201	gforget	1.3	if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
202			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
203	gforget	1.8	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
204	gforget	1.3	myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'};
205			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
206
207			if multiBasins;
208	gforget	1.8
209			%meridional streamfunction in Atlantic:
210			fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN; title0='Atlantic Meridional Stream Function';
211	gforget	1.3	if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
212			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
213	gforget	1.8	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
214	gforget	1.3	myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'};
215			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
216	gforget	1.8
217			%meridional streamfunction Pacific+Indian:
218			fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN; title0='Pac+Ind Meridional Stream Function';
219	gforget	1.3	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	gforget	1.8	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
222	gforget	1.3	myCaption={myYmeanTxt,'mean -- Pac+Ind overturning streamfunction (Sv)'};
223			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
224	gforget	1.8
225	gforget	1.3	end;
226
227	gforget	1.1
228			end;