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

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

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


Revision 1.4 - (hide annotations) (download)
Sat Dec 29 17:58:16 2012 UTC (12 years, 6 months ago) by gforget
Branch: MAIN
Changes since 1.3: +231 -3 lines
- add plotting sequence (ecco_v4/basic_diags_ecco_disp.m codes) in the
  various diagnostics sets, selected by the caller setting userStep==-1;

1 gforget 1.1
2     if userStep==1;%diags to be computed
3 gforget 1.2 %here _s stands for cumulative sum and _2s for cumulative sum of squares
4     %the _s should be stated first
5     listDiags=['THETA_s SALT_s SIatmQnt_s SIatmFW_s oceQnet_s oceFWflx_s ' ...
6     'fldTZ_s fldTM_s curlTau_s fldETAN_s fldETANLEADS_s fldMLD_s ' ...
7     'THETA_2s SALT_2s SIatmQnt_2s SIatmFW_2s oceQnet_2s oceFWflx_2s ' ...
8     'fldTZ_2s fldTM_2s curlTau_2s fldETAN_2s fldETANLEADS_2s fldMLD_2s'];
9 gforget 1.4
10 gforget 1.1 elseif userStep==2;%input files and variables
11     listFlds={ 'THETA','SALT','SIatmQnt','SIatmFW ','oceQnet ','oceFWflx','oceTAUX','oceTAUY','ETAN','sIceLoad','MXLDEPTH'};
12     listFldsNames=deblank(listFlds);
13     listFiles={'state_3d_set1','state_2d_set1','other_2d_set1'};
14     listSubdirs={[dirModel 'diags/OTHER/' ],[dirModel 'diags/STATE/' ]};
15 gforget 1.4
16     elseif userStep==3;%computation
17 gforget 1.1 %mask fields:
18     SIatmQnt=SIatmQnt.*mygrid.mskC(:,:,1);
19     SIatmFW=SIatmFW.*mygrid.mskC(:,:,1);
20     oceQnet=oceQnet.*mygrid.mskC(:,:,1);
21     oceFWflx=oceFWflx.*mygrid.mskC(:,:,1);
22     fldTX=oceTAUX.*mygrid.mskW(:,:,1);
23     fldTY=oceTAUY.*mygrid.mskS(:,:,1);
24     %compute Eastward/Northward wind stresses:
25     [fldTZ,fldTM]=calc_UEVNfromUXVY(fldTX,fldTY);
26     %compute wind stress curl:
27     curlTau=calc_UV_curl(fldTX, fldTY,1 );%the doMask argument should not matter as msk was already applied
28     %mask and re-arrange fields:
29     fldETAN=ETAN.*mygrid.mskC(:,:,1);
30     fldETANLEADS=(ETAN+sIceLoad/myparms.rhoconst).*mygrid.mskC(:,:,1);
31     fldMLD=MXLDEPTH.*mygrid.mskC(:,:,1);
32     %
33     THETA=THETA.*mygrid.mskC;
34     SALT=SALT.*mygrid.mskC;
35 gforget 1.2 %
36     if ii>1;
37     fileMatPrev=['diags_set_' tmp1 '_' num2str(listTimes(ii-1)) '.mat'];
38 gforget 1.4 listTimesBak=listTimes;
39 gforget 1.2 load([dirMat fileMatPrev]);
40 gforget 1.4 listTimes=listTimesBak;
41 gforget 1.2 end;
42     for jj=1:length(listDiags)/2;
43     myDiag=listDiags{jj}(1:end-2);
44     if ii==1;
45     eval([myDiag '_s=0*' myDiag ';']);
46     eval([myDiag '_2s=0*' myDiag '.^2;']);
47     end;
48     eval([myDiag '_s=' myDiag '_s+' myDiag ';']);
49     eval([myDiag '_2s=' myDiag '_2s+' myDiag '.^2;']);
50     end;
51    
52 gforget 1.4 %===================== COMPUTATIONAL SEQUENCE ENDS =========================%
53     %===================== PLOTTING SEQUENCE BEGINS =========================%
54    
55     elseif userStep==0;%loading / post-processing of mat files
56    
57     %load last cumsum
58     alldiag=alldiag_load(dirMat,[fileMat '_*.mat'],'',myparms.recInAve(2));
59     %load first cumsum
60     if myparms.recInAve(1)>1;
61     tmpdiag=alldiag_load(dirMat,[fileMat '_*.mat'],'',myparms.recInAve(1)-1);
62     for ii=1:length(alldiag.listDiags);
63     tmp0=alldiag.listDiags{ii};
64     if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps');
65     tmp1=getfield(alldiag,alldiag.listDiags{ii});
66     tmp2=getfield(tmpdiag,alldiag.listDiags{ii});
67     alldiag=setfield(alldiag,tmp0,tmp1-tmp2);
68     end;
69     end;
70     end;
71     %accomodate missing fields (old MITgcm version)
72     listFlds={'SIatmQnt_s','oceQnet_s','SIatmFW_s','oceFWflx_s',...
73     'fldTZ_s','fldTM_s','curlTau_s','fldETAN_s','fldETANLEADS_s',...
74     'SIatmQnt_2s','oceQnet_2s','SIatmFW_2s','oceFWflx_2s',...
75     'fldTZ_2s','fldTM_2s','curlTau_2s','fldETAN_2s','fldETANLEADS_2s'};
76     for ii=1:length(listFlds);
77     if ~sum(strcmp(fieldnames(alldiag),listFlds{ii}));
78     eval(['alldiag.' listFlds{ii} '=NaN*alldiag.fldETAN_s;']);
79     end;
80     end;
81     %
82     n=diff(myparms.recInAve)+1;
83     for ii=1:length(alldiag.listDiags)/2;
84     tmp0=alldiag.listDiags{ii};
85     if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps');
86     tmp1=1/n*getfield(alldiag,tmp0);
87     tmp2=1/n*getfield(alldiag,[tmp0(1:end-1) '2s']);
88     tmp2=(tmp2-tmp1.^2);
89     tmp2=n/(n-1)*tmp2;
90     %tmp2(tmp2<0)=0;
91     tmp2=sqrt(tmp2);
92     alldiag=setfield(alldiag,[tmp0(1:end-2) '_mean'],tmp1);
93     alldiag=setfield(alldiag,[tmp0(1:end-2) '_std'],tmp2);
94     end;
95     end;
96    
97     diagsWereLoaded=1
98    
99     elseif userStep==-1;%plotting
100    
101     if isempty(setDiagsParams);
102     choicePlot={'all'};
103     else;
104     choicePlot=setDiagsParams;
105     end;
106    
107     if sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'qnet'));
108    
109     %qnet from ocean+ice:
110     fld=-alldiag.SIatmQnt_mean;
111     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]; title0='QNET to ocean+ice';
112     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*10; end;
113     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
114     myCaption={myYmeanTxt,'mean -- QNET to ocean+ice (W/m2)'};
115     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
116    
117     %qnet to ocean:
118     fld=alldiag.oceQnet_mean;
119     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]; title0='QNET to ocean';
120     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*10; end;
121     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
122     myCaption={myYmeanTxt,'mean -- QNET to ocean (W/m2)'};
123     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
124    
125     if multiTimes;
126     %qnet from ocean+ice:
127     fld=alldiag.SIatmQnt_std;
128     cc=[[0:5:25] 35 [50:25:200] [250 300]]; title0='std(QNET to ocean+ice)';
129     if doAnomalies; cc=scaleAnom*[0:0.1:1]*5; end;
130     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
131     myCaption={myYmeanTxt,' standard deviation -- QNET to ocean+ice (W/m2)'};
132     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
133    
134     %qnet from ocean:
135     fld=alldiag.oceQnet_std;
136     cc=[[0:5:25] 35 [50:25:200] [250 300]]; title0='std(QNET to ocean)';
137     if doAnomalies; cc=scaleAnom*[0:0.1:1]*5; end;
138     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
139     myCaption={myYmeanTxt,' standard deviation -- QNET to ocean (W/m2)'};
140     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
141     end;
142    
143     end;
144    
145    
146     if sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'fwf'));
147    
148     %FW flux from ocean+ice:
149     fld=-alldiag.SIatmFW_mean/1000;%conversion to m/s
150     fld=fld*86400*1000;%conversion to mm/day
151     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]*0.06; title0='EMPMR to ocean+ice';
152     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.5; end;
153     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
154     myCaption={myYmeanTxt,'mean -- EMPMR from ocean+ice (mm/day)'};
155     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
156    
157     %FW flux from ocean:
158     fld=-alldiag.oceFWflx_mean/1000;%conversion to m/s
159     fld=fld*86400*1000;%conversion to mm/day
160     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]*0.06; title0='EMPMR to ocean';
161     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.5; end;
162     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
163     myCaption={myYmeanTxt,'mean -- EMPMR from ocean (mm/day)'};
164     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
165    
166     if multiTimes;
167     %empmr from ocean+ice:
168     fld=alldiag.SIatmFW_std*86400;%conversion to mm/day
169     cc=[[0:5:25] 35 [50:25:200] [250 300]]*0.04; title0='std(EMPMR to ocean+ice)';
170     if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.5; end;
171     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
172     myCaption={myYmeanTxt,' standard deviation -- EMPMR to ocean+ice (W/m2)'};
173     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
174    
175     %empmr from ocean:
176     fld=alldiag.oceFWflx_std*86400;%conversion to mm/day
177     cc=[[0:5:25] 35 [50:25:200] [250 300]]*0.04; title0='std(EMPMR to ocean)';
178     if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.5; end;
179     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
180     myCaption={myYmeanTxt,' standard deviation -- EMPMR to ocean (W/m2)'};
181     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
182     end;
183    
184     end;
185    
186     if sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'tau'));
187    
188     %zonal wind stress:
189     fld=alldiag.fldTZ_mean;
190     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='zonal wind stress';
191     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
192     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
193     myCaption={myYmeanTxt,'mean -- zonal wind stress (N/m2)'};
194     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
195    
196     %meridional wind stress:
197     fld=alldiag.fldTM_mean;
198     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='meridional wind stress';
199     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
200     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
201     myCaption={myYmeanTxt,'mean -- meridional wind stress (N/m2)'};
202     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
203    
204     fld=alldiag.curlTau_mean;
205     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/5e8; title0='wind stress curl';
206     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
207     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
208     myCaption={myYmeanTxt,'mean -- wind stress curl (N/m3)'};
209     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
210    
211     if multiTimes;
212     %zonal wind stress:
213     fld=alldiag.fldTZ_std;
214     cc=[[0:5:25] 35 [50:25:200] [250 300]]/2000; title0='std(tauZ)';
215     if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.005; end;
216     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
217     myCaption={myYmeanTxt,' standard deviation -- tauZ (W/m2)'};
218     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
219    
220     %meridional wind stress:
221     fld=alldiag.fldTM_std;
222     cc=[[0:5:25] 35 [50:25:200] [250 300]]/2000; title0='std(tauM)';
223     if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.005; end;
224     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
225     myCaption={myYmeanTxt,' standard deviation -- tauM (W/m2)'};
226     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
227    
228     fld=alldiag.curlTau_std;
229     cc=[[0:5:25] 35 [50:25:200] [250 300]]/1e9; title0='wind stress curl';
230     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.01; end;
231     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
232     myCaption={myYmeanTxt,'standard deviation -- tauCurl (N/m3)'};
233     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
234     end;
235    
236     end;
237    
238     if sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'ssh'));
239    
240     %ETAN:
241     fld=alldiag.fldETAN_mean;
242     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/100; title0='sea surface height (EXCLUDING ice)';
243     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.05; end;
244     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
245     myCaption={myYmeanTxt,'mean -- sea surface height (EXCLUDING ice, in m)'};
246     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
247    
248     %ETANLEADS:
249     fld=alldiag.fldETANLEADS_mean;
250     cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/100; title0='sea surface height (INCLUDING ice)';
251     if doAnomalies; cc=scaleAnom*[-1:0.1:1]*0.05; end;
252     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
253     myCaption={myYmeanTxt,'mean -- sea surface height (INCLUDING ice, in m)'};
254     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
255    
256     if multiTimes;
257     %ETAN:
258     fld=alldiag.fldETAN_std;
259     cc=[0:25:500]/2500; title0='std(ETAN)';
260     if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.02; end;
261     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
262     myCaption={myYmeanTxt,' standard deviation -- sea surface height (EXCLUDING ice, in m)'};
263     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
264    
265     %ETANLEADS:
266     fld=alldiag.fldETANLEADS_std;
267     cc=[0:25:500]/2500; title0='std(ETANLEADS)';
268     if doAnomalies; cc=scaleAnom*[0:0.1:1]*0.02; end;
269     figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
270     myCaption={myYmeanTxt,' standard deviation -- sea surface height (INCLUDING ice, in m)'};
271     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
272     end;
273    
274     end;
275    
276 gforget 1.1 end;
277 gforget 1.4
278    

  ViewVC Help
Powered by ViewVC 1.1.22