/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_C.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_C.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (hide annotations) (download)
Mon Jan 7 23:26:43 2013 UTC (12 years, 6 months ago) by gforget
Branch: MAIN
Changes since 1.2: +4 -2 lines
- fixes : backward compatibility, anomaly caxes, anomaly computation.

1 gforget 1.1
2     if userStep==1;%diags to be computed
3     listDiags='fldTzonmean fldSzonmean fldETANzonmean fldETANLEADSzonmean fldSIzonmean fldMLDzonmean';
4     listDiags=[listDiags ' IceAreaNorth IceAreaSouth IceVolNorth IceVolSouth SnowVolNorth SnowVolSouth'];
5     listDiags=[listDiags ' Ueq Teq Tglo Sglo TgloProf SgloProf'];
6     elseif userStep==2;%input files and variables
7     listFlds={ 'THETA','SALT','ETAN','SIarea','SIheff','SIhsnow','sIceLoad','MXLDEPTH','UVELMASS','VVELMASS'};
8     listFldsNames=deblank(listFlds);
9     listFiles={'state_2d_set1','other_2d_set1','state_3d_set1','trsp_3d_set1'};
10     listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/']};
11     elseif userStep==3;%computational part;
12     %mask and re-arrange fields:
13     fldT=THETA.*mygrid.mskC; fldS=SALT.*mygrid.mskC;
14     fldETAN=ETAN.*mygrid.mskC(:,:,1);
15     fldETANLEADS=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
16     fldSI=SIarea.*mygrid.mskC(:,:,1);
17     fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
18     ux=UVELMASS.*mygrid.mskW; vy=VVELMASS.*mygrid.mskS;
19     [ue,un]=calc_UEVNfromUXVY(ux,vy);
20    
21     %compute zonal means:
22     [fldTzonmean]=calc_zonmean_T(fldT);
23     [fldSzonmean]=calc_zonmean_T(fldS);
24     [fldETANzonmean]=calc_zonmean_T(fldETAN);
25     [fldETANLEADSzonmean]=calc_zonmean_T(fldETANLEADS);
26     [fldSIzonmean]=calc_zonmean_T(fldSI);
27     [fldMLDzonmean]=calc_zonmean_T(fldMLD);
28    
29     %compute hemispheric ice sums:
30     fld=SIarea.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorth=sum(fld);
31     fld=SIarea.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouth=sum(fld);
32     fld=SIheff.*mygrid.RAC.*(mygrid.YC>0); IceVolNorth=sum(fld);
33     fld=SIheff.*mygrid.RAC.*(mygrid.YC<0); IceVolSouth=sum(fld);
34     fld=SIhsnow.*mygrid.RAC.*(mygrid.YC>0); SnowVolNorth=sum(fld);
35     fld=SIhsnow.*mygrid.RAC.*(mygrid.YC<0); SnowVolSouth=sum(fld);
36    
37     %equatorial sections of T and U:
38     [secX,secY,Teq]=gcmfaces_section([],0,fldT,1);
39     [secX,secY,Ueq]=gcmfaces_section([],0,ue,1);
40    
41     %global means and profiles:
42     msk=mygrid.mskC.*mygrid.hFacC;
43     msk=msk.*mk3D(mygrid.RAC,msk).*mk3D(mygrid.DRF,msk);
44     Tglo=nansum(fldT.*msk)./nansum(msk);
45     Sglo=nansum(fldS.*msk)./nansum(msk);
46     nr=length(mygrid.RC);
47     TgloProf=zeros(nr,1); SgloProf=zeros(nr,1);
48     for kk=1:nr;
49     rac=mygrid.RAC.*mygrid.mskC(:,:,kk);
50     TgloProf(kk)=nansum(rac.*fldT(:,:,kk))/nansum(rac);
51     SgloProf(kk)=nansum(rac.*fldS(:,:,kk))/nansum(rac);
52     end;
53 gforget 1.2
54    
55     %===================== COMPUTATIONAL SEQUENCE ENDS =========================%
56     %===================== PLOTTING SEQUENCE BEGINS =========================%
57    
58     elseif userStep==-1;%plotting
59    
60     if isempty(setDiagsParams);
61     choicePlot={'all'};
62     else;
63     choicePlot=setDiagsParams;
64     end;
65    
66     if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'zmtend')));
67    
68     figureL; set(gcf,'Renderer','zbuffer');
69     %cc=[-2:0.25:2];
70     cc=1/100*[[-200:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:200]]
71     if doAnomalies; cc=scaleAnom/10*cc; end;
72     X=mygrid.LATS*ones(1,length(mygrid.RC)); Y=ones(length(mygrid.LATS),1)*(mygrid.RC');
73     %T
74     subplot(2,1,1);
75     fld=annualmean(alldiag.listTimes,alldiag.fldTzonmean,'last')-annualmean(alldiag.listTimes,alldiag.fldTzonmean,'first');
76     depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
77     cbar=gcmfaces_cmap_cbar(cc); title('zonal mean T anomaly');
78     %S
79     subplot(2,1,2);
80     fld=annualmean(alldiag.listTimes,alldiag.fldSzonmean,'last')-annualmean(alldiag.listTimes,alldiag.fldSzonmean,'first');
81     depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
82     cbar=gcmfaces_cmap_cbar(cc/4); title('zonal mean S anomaly');
83     %
84     myCaption={myYmeanTxt,', last year minus first year -- zonal mean temperature (degC; top) and salinity (psu; bottom)'};
85     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
86    
87     end;
88    
89     if sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'eq'));
90    
91     %time mean equator sections :
92     figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2];
93     [secX,secY,LONeq]=gcmfaces_section([],0,mygrid.XC,1);
94     X=LONeq*ones(1,length(mygrid.RC)); Y=ones(length(LONeq),1)*mygrid.RC';
95     X=circshift(X,[-180 0]); X(X<0)=X(X<0)+360;
96     %T
97     subplot(2,1,1);
98     fld=annualmean(alldiag.listTimes,alldiag.Teq,'all'); fld=circshift(fld,[-180 0]);
99     depthStretchPlot('pcolor',{X,Y,fld},[0:50:350],[0 350 350]); shading interp;
100     if ~doAnomalies; cc=18+6*cc0; else; cc=scaleAnom/10*cc0*2; end;
101     cbar=gcmfaces_cmap_cbar(cc); title('equator T (degC)');
102     %U
103     subplot(2,1,2);
104     fld=annualmean(alldiag.listTimes,alldiag.Ueq,'all'); fld=circshift(fld,[-180 0]);
105     depthStretchPlot('pcolor',{X,Y,fld},[0:50:350],[0 350 350]); shading interp;
106     if ~doAnomalies; cc=cc0/5*2; else; cc=scaleAnom/10*cc0/5; end;
107     cbar=gcmfaces_cmap_cbar(cc); title('equator zonal velocity (m/s)');
108     %
109     myCaption={myYmeanTxt,' mean -- equator temperature (degC;top) and zonal velocity (m/s;bottom)'};
110     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
111    
112     end;
113    
114     if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'gmtime')));
115    
116     %global mean T/S:
117     figureL;
118     subplot(2,1,1); vec=alldiag.Tglo;
119     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
120     aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
121     legend('monthly','annual mean'); title('Global Mean Temperature (degC)');
122     subplot(2,1,2); vec=alldiag.Sglo;
123     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2);
124     aa=axis; axis([min(TT) max(TT) aa(3:4)]); grid on;
125     legend('monthly','annual mean'); title('Global Mean Salinity (psu)');
126     myCaption={'global mean T (degC; top) and S (psu; bottom)'};
127     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
128    
129     %global mean T/S profiles:
130     figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2];
131     X=TT*ones(1,length(mygrid.RC)); Y=ones(length(TT),1)*(mygrid.RC'); X=X'; Y=Y';
132     %T
133     subplot(2,1,1);
134     [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.TgloProf),'first');
135     depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
136     if ~doAnomalies; cc=cc0/5; else; cc=scaleAnom/10*cc0/10; end;
137     cbar=gcmfaces_cmap_cbar(cc); title('global mean T minus first year (K)');
138     %S
139     subplot(2,1,2);
140     [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.SgloProf),'first');
141     depthStretchPlot('pcolor',{X,Y,fld}); shading interp;
142     if ~doAnomalies; cc=cc0/50; else; cc=scaleAnom/10*cc0/50; end;
143     cbar=gcmfaces_cmap_cbar(cc); title('global mean S minus first year (psu)');
144     %
145     myCaption={'global mean temperature (K; top) and salinity (psu; bottom) minus first year'};
146     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
147    
148     end;
149    
150    
151     if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'lzmtime')));
152    
153     listLatPlot=[-65 -25 0 25 65]
154     %zonal mean temperature profile at selected latitude:
155     for iLatPlot=1:length(listLatPlot);
156     tmp1=abs(mygrid.LATS-listLatPlot(iLatPlot));
157     iLat=find(tmp1==min(tmp1)); iLat=iLat(1);
158     figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2];
159     if doAnomalies; cc=scaleAnom/10*cc; end;
160     %T
161     subplot(2,1,1);
162     [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.fldTzonmean(iLat,:,:)),'first');
163     X=TT*ones(1,length(mygrid.RC)); Y=ones(length(TT),1)*(mygrid.RC'); X=X'; Y=Y';
164     title0=['mean T minus first year (K) at lat ~ ' num2str(listLatPlot(iLatPlot))];
165     depthStretchPlot('pcolor',{X,Y,fld}); shading interp; cbar=gcmfaces_cmap_cbar(cc/5); title(title0);
166     %S
167     subplot(2,1,2);
168     [tmp1,fld]=annualmean(alldiag.listTimes,squeeze(alldiag.fldSzonmean(iLat,:,:)),'first');
169     X=TT*ones(1,length(mygrid.RC)); Y=ones(length(TT),1)*(mygrid.RC'); X=X'; Y=Y';
170     title0=['mean S minus first year (psu) at lat ~ ' num2str(listLatPlot(iLatPlot))];
171     depthStretchPlot('pcolor',{X,Y,fld}); shading interp; cbar=gcmfaces_cmap_cbar(cc/25); title(title0);
172    
173     myCaption={'mean temperature (top; K) and salinity (bottom; psu) minus first year ',...
174     ['at lat $\\approx$ ' num2str(listLatPlot(iLatPlot))]};
175     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
176     end;
177    
178     end;
179    
180    
181     if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'szmtime')));
182    
183     %zonal mean SST/SSS:
184     figureL; set(gcf,'Renderer','zbuffer'); cc0=[-2:0.25:2]; kk=1;
185     x=TT*ones(1,length(mygrid.LATS)); y=ones(nt,1)*mygrid.LATS'; x=x'; y=y';
186     %T
187     subplot(2,1,1);
188     fld=squeeze(alldiag.fldTzonmean(:,kk,:)); [tmp1,fld]=annualmean(alldiag.listTimes,fld,'first');
189     if ~doAnomalies; cc=cc0*3; else; cc=scaleAnom/10*cc0; end;
190     pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]); cbar=gcmfaces_cmap_cbar(cc);
191     title('zonal mean SST minus first year (degC)');
192     %S
193     subplot(2,1,2);
194     fld=squeeze(alldiag.fldSzonmean(:,kk,:)); [tmp1,fld]=annualmean(alldiag.listTimes,fld,'first');
195     if ~doAnomalies; cc=cc0*2/5; else; cc=scaleAnom/10*cc0/2; end;
196     pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]); cbar=gcmfaces_cmap_cbar(cc);
197     title('zonal mean SSS minus first year (psu)');
198     %
199     myCaption={'zonal mean temperature (degC; top) and salinity ',...
200     '(psu; bottom) minus first year (psu)',[' at ' num2str(-mygrid.RC(kk)) 'm depth']};
201     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
202    
203     figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2]; kk=1;
204     if doAnomalies; cc=scaleAnom/10*cc/2; end;
205     x=TT*ones(1,length(mygrid.LATS)); y=ones(nt,1)*mygrid.LATS'; x=x'; y=y';
206     %SSH
207     subplot(2,1,1);
208     [tmp1,fld]=annualmean(alldiag.listTimes,alldiag.fldETANLEADSzonmean,'first');
209     pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
210     cbar=gcmfaces_cmap_cbar(2*cc/25); title(['zonal mean SSH minus first year (INCLUDING ice, in m)']);
211     %ETAN
212     subplot(2,1,2);
213     [tmp1,fld]=annualmean(alldiag.listTimes,alldiag.fldETANzonmean,'first');
214     pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
215     cbar=gcmfaces_cmap_cbar(4*cc/5); title(['zonal mean SSH minus first year (EXCLUDING ice, in m)']);
216     %
217     myCaption={'zonal mean SSH (m) minus first year, including ice (top) and below ice (bottom) '};
218     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
219    
220     figureL; set(gcf,'Renderer','zbuffer'); cc=[-2:0.25:2]; kk=1;
221     if doAnomalies; cc=scaleAnom/10*cc; end;
222     x=TT*ones(1,length(mygrid.LATS)); y=ones(nt,1)*mygrid.LATS'; x=x'; y=y';
223     %SSH
224     subplot(2,1,1);
225     fld=squeeze(alldiag.fldSIzonmean(:,tt));
226     pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
227 gforget 1.3 if doAnomalies; cc2=scaleAnom/50*cc; else; cc2=0.5+cc/4; end;
228     cbar=gcmfaces_cmap_cbar(cc2); title(['zonal mean Ice conc. -- in m ']);
229 gforget 1.2 %MLD
230     subplot(2,1,2);
231     fld=squeeze(alldiag.fldMLDzonmean(:,tt));
232     pcolor(x,y,fld); shading interp; axis([TT(1) TT(end) -90 90]);
233 gforget 1.3 if doAnomalies; cc2=scaleAnom*40*cc; else; cc2=200+500*cc/5; end;
234     cbar=gcmfaces_cmap_cbar(cc2); title(['zonal mean MLD -- in m ']);
235 gforget 1.2 %
236     myCaption={'zonal mean ice concentration (no units) and mixed layer depth (m)'};
237     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
238    
239     end;
240    
241    
242     if multiTimes&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'icetime')));
243    
244     %sea ice area
245     figureL;
246     subplot(2,1,1); vec=alldiag.IceAreaNorth/1e12;
247     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 20]);
248     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
249     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice cover (in 10^12m^2)');
250     subplot(2,1,2); vec=alldiag.IceAreaSouth/1e12;
251     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 25]);
252     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice cover (in 10^12m^2)');
253     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
254     myCaption={'sea ice cover (in $10^12m^2$) in northern (top) and southern (bottom) hemisphere'};
255     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
256    
257     %sea ice volume
258     figureL;
259     subplot(2,1,1); vec=alldiag.IceVolNorth/1e12;
260     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 50]);
261     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice volume (in 10^12m^3)');
262     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]*2; axis(aa); end;
263     subplot(2,1,2); vec=alldiag.IceVolSouth/1e12;
264     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 15]);
265     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice volume (in 10^12m^3)');
266     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]*2; axis(aa); end;
267     myCaption={'sea ice volume (in $10^12m^3$) in northern (top) and southern (bottom) hemisphere'};
268     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
269    
270     %snow volume
271     figureL;
272     subplot(2,1,1); vec=alldiag.SnowVolNorth/1e12;
273     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 5]);
274     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere snow volume (in 10^12m^3)');
275     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
276     subplot(2,1,2); vec=alldiag.SnowVolSouth/1e12;
277     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 4]);
278     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere snow volume (in 10^12m^3)');
279     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
280     myCaption={'snow volume (in $10^12m^3$) in northern (top) and southern (bottom) hemisphere'};
281     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
282    
283     if ~doAnomalies; %does not work for anomalies -> need to compute ratio in basic_diags_ecco
284     %sea ice thickness
285     figureL;
286     subplot(2,1,1); vec=alldiag.IceVolNorth./alldiag.IceAreaNorth;
287     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 5]);
288     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere ice thickness (in m)');
289     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
290     subplot(2,1,2); vec=alldiag.IceVolSouth./alldiag.IceAreaSouth;
291     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 2]);
292     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere ice thickness (in m)');
293     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
294     myCaption={'sea ice thickness (in m) in northern (top) and southern (bottom) hemisphere'};
295     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
296    
297     %sea ice thickness
298     figureL;
299     subplot(2,1,1); vec=alldiag.SnowVolNorth./alldiag.IceAreaNorth;
300     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 0.4]);
301     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Northern Hemisphere snow thickness (in m)');
302     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
303     subplot(2,1,2); vec=alldiag.SnowVolSouth./alldiag.IceAreaSouth;
304     plot(TT,vec,'LineWidth',2); hold on; plot(TT,runmean(vec,myNmean,2),'r','LineWidth',2); axis([min(TT) max(TT) 0 0.25]);
305     grid on; if myNmean>0; legend('monthly','annual mean'); end; title('Southern Hemisphere snow thickness (in m)');
306     if doAnomalies; aa(3:4)=scaleAnom*[-1 1]; axis(aa); end;
307     myCaption={'snow thickness (in m) in northern (top) and southern (bottom) hemisphere'};
308     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
309     end;
310    
311     end;
312    
313 gforget 1.1 end;

  ViewVC Help
Powered by ViewVC 1.1.22