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
2	layersOffline=1;
3	layersEulerian=0;
4	integrateFromBottom=1;
5	if isempty(whos('setDiagsParams')); setDiagsParams={}; end;
6
7	if length(setDiagsParams)==2;
8	layersName=setDiagsParams{1};
9	layersLims=setDiagsParams{2};
10	layersOffline=1;
11	layersEulerian=1;%hacked : should be 0
12	elseif length(setDiagsParams)==1;
13	layersName=setDiagsParams{1};
14	else;
15	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	end;
21
22	layersFileSuff=['LAYERS_' layersName];
23	if layersOffline; layersFileSuff=['LAYERS_offline_' layersName]; end;
24
25	%override default file name:
26	%---------------------------
27	fileMat=['diags_set_' layersFileSuff];
28
29	if userStep==1;%diags to be computed
30	listDiags=['layersParams gloOV gloThick gloDelThick gloBAR gloMT gloD'];
31	listBasins=1;
32	if sum([90 1170]~=mygrid.ioSize)==0;
33	listBasins=[1:3];
34	mskC=v4_basin({'atlExt'});
35	if isempty(mskC); listBasins=1; end;
36	end;
37	if length(listBasins)==3;
38	listDiags=[listDiags ' atlOV atlThick atlDelThick atlBAR atlMT atlD'];
39	listDiags=[listDiags ' pacindOV pacindThick pacindDelThick pacindBAR pacindMT pacindD'];
40	end;
41	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	listSubdirs={[dirModel 'diags/LAYERS/' ],[dirModel 'diags/']};
47	%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	elseif userStep==3;
54
55	%override default file name:
56	%---------------------------
57	fileMat=['diags_set_' layersFileSuff '_' num2str(tt) '.mat'];
58
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	fldD=[];
68	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
81	if strcmp(layersName,'theta');
82	tracer=THETA;
83	if isempty(whos('layersLims')); layersLims=[-2:35]; end;
84	elseif strcmp(layersName,'salt');
85	tracer=SALT;
86	if isempty(whos('layersLims')); layersLims=[33:0.1:38]; end;
87	elseif strcmp(layersName,'sigma');
88	P=mk3D(-mygrid.RC,THETA);
89	t=convert2vector(THETA);
90	s=convert2vector(SALT);
91	p=convert2vector(P);
92	pref=pref0+0*p;
93	[rhop,rhpis,rhor] = density(t(:),s(:),p(:),pref(:));
94	rhor=reshape(rhor,size(t));
95	tracer=convert2vector(rhor)-1000;
96	if isempty(whos('layersLims')); layersLims=layersLims0; end;
97	end
98
99	layersGrid=(layersLims(1:end-1)+layersLims(2:end))'/2;
100	[fldU,fldV]=layers_remap({U,V},'extensive',tracer,layersGrid,2);
101	fldHw=[]; fldHs=[];
102	[fldD]=layers_remap(P,'intensive',tracer,layersGrid,2);
103
104	end;%if ~layersOffline;
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	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
116	%the overturning streamfunction computation itself
117	layersOV=calc_overturn(fldU.mskW3d,fldV.mskS3d,1,{'dh'});
118
119	%the associated barotropic streamfunction (for checking)
120	layersBAR=calc_barostream(fldU.mskW3d,fldV.mskS3d,1,{'dh'});
121
122	%meridional transport per layer:
123	layersMT=diff(layersOV,1,2);
124
125	if ~isempty(fldD);
126	%compute zonal mean depth from fldD
127	layersD=calc_zonmean_T(fldD.*mskC3d);
128	else;
129	layersD=[];
130	end;
131
132	%the associated thickness
133	if isempty(fldHw);
134	layersThick=[];
135	layersDelThick=[];
136	else;
137	%interpolate to tracer points
138	fldH=0*fldHw;
139	[fldHwp,fldHsp]=exch_UV(LaHw1SLT.mskW3d,LaHs1SLT.mskS3d);
140	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
163	%store to global, atlantic or Pac+Ind arrays:
164	if bb==1;
165	gloOV=layersOV; gloThick=layersThick; gloDelThick=layersDelThick;
166	gloBAR=layersBAR; gloMT=layersMT; gloD=layersD;
167	elseif bb==2;
168	kk=find(mygrid.LATS<-35\|mygrid.LATS>70);
169	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	elseif bb==3;
174	kk=find(mygrid.LATS<-35\|mygrid.LATS>65);
175	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	end;
180	end;
181
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
189	%===================== COMPUTATIONAL SEQUENCE ENDS =========================%
190	%===================== PLOTTING SEQUENCE BEGINS =========================%
191
192	elseif userStep==-1;%plotting
193
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
200	%meridional streamfunction (Eulerian) :
201	fld=mean(alldiag.gloOV(:,:,tt),3); fld(fld==0)=NaN; title0='Meridional Stream Function';
202	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	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
205	myCaption={myYmeanTxt,'mean -- overturning streamfunction (Sv)'};
206	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
207
208	if multiBasins;
209
210	%meridional streamfunction in Atlantic:
211	fld=mean(alldiag.atlOV(:,:,tt),3); fld(fld==0)=NaN; title0='Atlantic Meridional Stream Function';
212	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	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
215	myCaption={myYmeanTxt,'mean -- Atlantic overturning streamfunction (Sv)'};
216	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
217
218	%meridional streamfunction Pacific+Indian:
219	fld=mean(alldiag.pacindOV(:,:,tt),3); fld(fld==0)=NaN; title0='Pac+Ind Meridional Stream Function';
220	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	pcolor(X,Y,fld); shading interp; cbar=gcmfaces_cmap_cbar(cc); title(title0); set(gca,'YDir','reverse');
223	myCaption={myYmeanTxt,'mean -- Pac+Ind overturning streamfunction (Sv)'};
224	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
225
226	end;
227
228	end;