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);
            %if isempty(whos('layersLims')); layersLims=1000+[30:0.1:38]'; end;
            if isempty(whos('layersLims')); layersLims=1000+[19:0.2:31]'; end;
        end
        
        layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
        [fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,2);
        fldHw=[]; fldHs=[];
        if 0;
          [fldD]=layers_remap(P,'intensive',tracer,layersGrid,2);
        end;

    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'});

    if 0;
    %meridional transport per layer:
    dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
    trU=fldU.*dyg; trV=fldV.*dxg;
    layersMT=zeros(length(mygrid.LATS_MASKS),length(layersGrid));   
      for iy=1:length(mygrid.LATS_MASKS);
        %get list ofpoints that form a zonal band:
        mskW=mygrid.LATS_MASKS(iy).mskWedge;
        vecW=gcmfaces_subset(mskW,trU);
        mskS=mygrid.LATS_MASKS(iy).mskSedge;
        vecS=gcmfaces_subset(mskS,trV);
        %store vertically integrated transport:
        layersMT(iy,:)=1e-6*(nansum(vecW,1)+nansum(vecS,1));
      end;
    else;
      layersMT=NaN*layersOV; 
    end;

    if 0;
      %...compute zonal mean depth from fldD
    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; gloDepth=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			tracer=convert2vector(rhor);
86	gforget	1.5	%if isempty(whos('layersLims')); layersLims=1000+[30:0.1:38]'; end;
87			if isempty(whos('layersLims')); layersLims=1000+[19:0.2:31]'; end;
88	gforget	1.1	end
89
90			layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
91	gforget	1.5	[fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,2);
92	gforget	1.1	fldHw=[]; fldHs=[];
93	gforget	1.5	if 0;
94			[fldD]=layers_remap(P,'intensive',tracer,layersGrid,2);
95			end;
96
97	gforget	1.1	end;%if ~layersOffline;
98	gforget	1.3
99			if sum([90 1170]~=mygrid.ioSize)==0;
100			listBasins=[1:3];
101			else;
102			listBasins=1;
103			end;
104
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			if 0;
122			%meridional transport per layer:
123			dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
124			trU=fldU.dyg; trV=fldV.dxg;
125			layersMT=zeros(length(mygrid.LATS_MASKS),length(layersGrid));
126			for iy=1:length(mygrid.LATS_MASKS);
127			%get list ofpoints that form a zonal band:
128			mskW=mygrid.LATS_MASKS(iy).mskWedge;
129			vecW=gcmfaces_subset(mskW,trU);
130			mskS=mygrid.LATS_MASKS(iy).mskSedge;
131			vecS=gcmfaces_subset(mskS,trV);
132			%store vertically integrated transport:
133			layersMT(iy,:)=1e-6*(nansum(vecW,1)+nansum(vecS,1));
134			end;
135			else;
136			layersMT=NaN*layersOV;
137			end;
138
139			if 0;
140			%...compute zonal mean depth from fldD
141			else;
142			layersD=[];
143			end;
144
145	gforget	1.1	%the associated thickness
146			if isempty(fldHw);
147			layersThick=[];
148			layersDelThick=[];
149			else;
150			%interpolate to tracer points
151			fldH=0*fldHw;
152	gforget	1.3	[fldHwp,fldHsp]=exch_UV(LaHw1SLT.mskW3d,LaHs1SLT.mskS3d);
153	gforget	1.1	for iF=1:fldH.nFaces;
154			tmpA=fldHwp{iF}(2:end,:,:);
155			tmpB=fldHwp{iF}(1:end-1,:,:);
156			tmpC=fldHsp{iF}(:,2:end,:);
157			tmpD=fldHsp{iF}(:,1:end-1,:);
158			tmpTot=tmpA+tmpB+tmpC+tmpD;
159			tmpNb=1(tmpA~=0)+1(tmpB~=0)+1(tmpC~=0)+1(tmpD~=0);
160			jj=find(tmpNb>0); tmpTot(jj)=tmpTot(jj)./tmpNb(jj);
161			jj=find(isnan(tmpTot)); tmpTot(jj)=0;
162			fldH{iF}=tmpTot;
163			end;
164			%compute zonal mean
165			fldH=calc_zonmean_T(fldH.*mskC3d);
166			%integrate to overturning points
167			if integrateFromBottom;
168			layersThick=[flipdim(cumsum(flipdim(fldH,2),2),2) zeros(size(fldH,1),1)];
169			else;
170			layersThick=[zeros(size(fldH,1),1) cumsum(fldH,2)];
171			end;
172			%compute dH/dLayer
173			layersDelThick=[zeros(size(fldH,1),1) fldH./( ones(size(fldH,1),1)*diff(layersLims') )];
174			end;
175	gforget	1.3
176			%store to global, atlantic or Pac+Ind arrays:
177			if bb==1;
178	gforget	1.5	gloOV=layersOV; gloThick=layersThick; gloDelThick=layersDelThick;
179			gloBAR=layersBAR; gloMT=layersMT; gloDepth=layersD;
180	gforget	1.3	elseif bb==2;
181			kk=find(mygrid.LATS<-35\|mygrid.LATS>70);
182	gforget	1.5	layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
183			if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
184			atlOV=layersOV; atlThick=layersThick; atlDelThick=layersDelThick;
185			atlBAR=layersBAR; atlMT=layersMT; atlD=layersD;
186	gforget	1.3	elseif bb==3;
187			kk=find(mygrid.LATS<-35\|mygrid.LATS>65);
188	gforget	1.5	layersOV(kk,:)=NaN; layersMT(kk,:)=NaN;
189			if ~isempty(layersThick); layersThick(kk,:)=NaN; layersDelThick(kk,:)=NaN; end;
190			pacindOV=layersOV; pacindThick=layersThick; pacindDelThick=layersDelThick;
191			pacindBAR=layersBAR; pacindMT=layersMT; pacindD=layersD;
192	gforget	1.3	end;
193			end;
194	gforget	1.1
195			layersParams.name=layersName;
196			layersParams.LC=layersGrid;
197			layersParams.LF=layersLims;
198			layersParams.isOffline=layersOffline;
199			layersParams.suffix=layersFileSuff;
200			layersParams.isEulerian=layersEulerian;
201	gforget	1.3
202
203			%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
204			%===================== PLOTTING SEQUENCE BEGINS =========================%
205
206			elseif userStep==-1;%plotting
207
208
209			%meridional streamfunction (Eulerian) :
210			fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN;
211			X=mygrid.LATSones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)(alldiag.layersParams.LF');
212			cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Meridional Stream Function';
213			if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
214			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
215			pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0);
216			myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'};
217			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
218
219			if multiBasins;
220			%meridional streamfunction (Eulerian):
221			fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN;
222			X=mygrid.LATSones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)(alldiag.layersParams.LF');
223			cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Atlantic Meridional Stream Function';
224			if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
225			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
226			pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0);
227			myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'};
228			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
229			%meridional streamfunction (residual):
230			fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN;
231			X=mygrid.LATSones(1,length(alldiag.layersParams.LF)); Y=ones(length(mygrid.LATS),1)(alldiag.layersParams.LF');
232			cc=[[-50:10:-30] [-24:3:24] [30:10:50]]; title0='Pac+Ind Meridional Stream Function';
233			if doAnomalies; cc=scaleAnom*[-1:0.1:1]; end;
234			figureL; set(gcf,'Renderer','zbuffer'); %set(gcf,'Units','Normalized','Position',[0.05 0.1 0.4 0.8]);
235			pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0);
236			myCaption={myYmeanTxt,'mean -- Pac+Ind overturning streamfunction (Sv)'};
237			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
238			end;
239
240	gforget	1.1
241			end;