Subversion Repositories pvinversion.ecmwf

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
% -------------------------------------------------------------------------
2
% Load files
3
% -------------------------------------------------------------------------
4
 
5
% Base directory and filename
6
base = '/net/rossby/lhome/sprenger/PV_Inversion_Tool/run/';
7
folder   = base;
8
filename = 'ANO_20060115_18';
9
disp([folder filename])
10
 
11
% First image (otherwise first image is not correctly written)
12
figname = [base '/test.eps'];
13
close;
14
fh=figure('Units','pixels','Position',[100 100 900 900])
15
set(gcf, 'PaperPosition', [2 1 15 10]);
16
print('-depsc2','-r0',figname);
17
 
18
% Load variables from file (on model levels)
19
z_th   = cdf_loadV(folder,filename,'TH');
20
z_uu   = cdf_loadV(folder,filename,'U');
21
z_vv   = cdf_loadV(folder,filename,'V');
22
z_qgpv = cdf_loadV(folder,filename,'QGPV');
23
z_psi  = cdf_loadV(folder,filename,'PSI');
24
 
25
% -------------------------------------------------------------------------
26
% Plot stream function
27
% -------------------------------------------------------------------------
28
 
29
for ilev=1:z_psi.nz
30
 
31
% Create a new figure
32
close;
33
fh=figure('Units','pixels','Position',[100 100 900 900])
34
 
35
% Set the geographical projection
36
m_proj('Equidistant Cylindrical','long',[z_psi.lonmin z_psi.lonmax],'lat',[z_psi.latmin z_psi.latmax]);
37
 
38
% Scale the plotting field for color map
39
fld=squeeze(z_psi.var(ilev,:,:));
40
c_map = scale_col(-4000:500:4000,fld);
41
 
42
% Plot Stream Function PSI
43
lat=z_psi.ymin + (0:z_psi.ny-1) * z_psi.dy;
44
lon=z_psi.xmin + (0:z_psi.nx-1) * z_psi.dx;
45
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
46
for icnt = 1: length(h)
47
  set( h(icnt), 'EdgeColor', 'none' )
48
end
49
 
50
% Add color bar
51
colormap('default');
52
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
53
caxis(c_map.caxis);
54
q=colorbar('vert');
55
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);
56
 
57
% Add the grid and the coast lines to the plot
58
m_grid;
59
m_coast('linewidth',1,'color','k');
60
title(num2str(z_psi.aklay(ilev)));
61
 
62
% Save figure
63
pre='';
64
if ( z_psi.aklay(ilev) < 10 )
65
  pre='0';
66
end
67
if ( z_psi.aklay(ilev) < 100 )
68
  pre=[pre '0' ];
69
end
70
if ( z_psi.aklay(ilev) < 1000 )
71
  pre=[pre '0' ];
72
end
73
if ( z_psi.aklay(ilev) < 10000 )
74
  pre=[pre '0' ];
75
end
76
 
77
figname = [ base '/sf_2d_' filename '_' pre num2str(z_psi.aklay(ilev)) '.eps' ];
78
set(gcf, 'PaperPosition', [2 1 15 10]);
79
print('-depsc2','-r0',figname);
80
 
81
% End loop over all levels
82
end
83
 
84
% -------------------------------------------------------------------------
85
% Plot anomaly of potential temperature
86
% -------------------------------------------------------------------------
87
 
88
for ilev=1:z_th.nz
89
 
90
% Create a new figure
91
close;
92
fh=figure('Units','pixels','Position',[100 100 900 900])
93
 
94
% Set the geographical projection
95
m_proj('Equidistant Cylindrical','long',[z_th.lonmin z_th.lonmax],'lat',[z_th.latmin z_th.latmax]);
96
 
97
% Scale the plotting field for color map
98
fld=squeeze(z_th.var(ilev,:,:));
99
c_map = scale_col(-20:2.5:20,fld);
100
 
101
% Plot potential temperature
102
lat=z_th.ymin + (0:z_th.ny-1) * z_th.dy;
103
lon=z_th.xmin + (0:z_th.nx-1) * z_th.dx;
104
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
105
for icnt = 1: length(h)
106
  set( h(icnt), 'EdgeColor', 'none' )
107
end
108
 
109
% Add color bar
110
colormap('default');
111
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
112
caxis(c_map.caxis);
113
q=colorbar('vert');
114
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);
115
 
116
% Add the grid and the coast lines to the plot
117
m_grid;
118
m_coast('linewidth',1,'color','k');
119
title(num2str(z_th.aklay(ilev)));
120
 
121
% Save figure
122
pre='';
123
if ( z_th.aklay(ilev) < 10 )
124
  pre='0';
125
end
126
if ( z_th.aklay(ilev) < 100 )
127
  pre=[pre '0' ];
128
end
129
if ( z_th.aklay(ilev) < 1000 )
130
  pre=[pre '0' ];
131
end
132
if ( z_th.aklay(ilev) < 10000 )
133
  pre=[pre '0' ];
134
end
135
 
136
figname = [ base '/ta_2d_' filename '_' pre num2str(z_th.aklay(ilev)) '.eps' ];
137
set(gcf, 'PaperPosition', [2 1 15 10]);
138
print('-depsc2','-r0',figname);
139
 
140
% End loop over all levels
141
end
142
 
143
% -------------------------------------------------------------------------
144
% Plot anomaly of meridional wind
145
% -------------------------------------------------------------------------
146
 
147
for ilev=1:z_vv.nz
148
 
149
% Create a new figure
150
close;
151
fh=figure('Units','pixels','Position',[100 100 900 900])
152
 
153
% Set the geographical projection
154
m_proj('Equidistant Cylindrical','long',[z_vv.lonmin z_vv.lonmax],'lat',[z_vv.latmin z_vv.latmax]);
155
 
156
% Scale the plotting field for color map
157
fld=squeeze(z_vv.var(ilev,:,:));
158
c_map = scale_col(-30:3:30,fld);
159
 
160
% Plot potential temperature
161
lat=z_vv.ymin + (0:z_vv.ny-1) * z_vv.dy;
162
lon=z_vv.xmin + (0:z_vv.nx-1) * z_vv.dx;
163
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
164
for icnt = 1: length(h)
165
  set( h(icnt), 'EdgeColor', 'none' )
166
end
167
 
168
% Add color bar
169
colormap('default');
170
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
171
caxis(c_map.caxis);
172
q=colorbar('vert');
173
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);
174
 
175
% Add the grid and the coast lines to the plot
176
m_grid;
177
m_coast('linewidth',1,'color','k');
178
title(num2str(z_vv.aklay(ilev)));
179
 
180
% Save figure
181
pre='';
182
if ( z_vv.aklay(ilev) < 10 )
183
  pre='0';
184
end
185
if ( z_vv.aklay(ilev) < 100 )
186
  pre=[pre '0' ];
187
end
188
if ( z_vv.aklay(ilev) < 1000 )
189
  pre=[pre '0' ];
190
end
191
if ( z_vv.aklay(ilev) < 10000 )
192
  pre=[pre '0' ];
193
end
194
 
195
figname = [ base '/va_2d_' filename '_' pre num2str(z_vv.aklay(ilev)) '.eps' ];
196
set(gcf, 'PaperPosition', [2 1 15 10]);
197
print('-depsc2','-r0',figname);
198
 
199
% End loop over all levels
200
end
201
 
202
% -------------------------------------------------------------------------
203
% Plot anomaly of zonal wind
204
% -------------------------------------------------------------------------
205
 
206
for ilev=1:z_uu.nz
207
 
208
% Create a new figure
209
close;
210
fh=figure('Units','pixels','Position',[100 100 900 900])
211
 
212
% Set the geographical projection
213
m_proj('Equidistant Cylindrical','long',[z_uu.lonmin z_uu.lonmax],'lat',[z_uu.latmin z_uu.latmax]);
214
 
215
% Scale the plotting field for color map
216
fld=squeeze(z_uu.var(ilev,:,:));
217
c_map = scale_col(-30:3:30,fld);
218
 
219
% Plot potential temperature
220
lat=z_uu.ymin + (0:z_uu.ny-1) * z_uu.dy;
221
lon=z_uu.xmin + (0:z_uu.nx-1) * z_uu.dx;
222
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
223
for icnt = 1: length(h)
224
  set( h(icnt), 'EdgeColor', 'none' )
225
end
226
 
227
% Add color bar
228
colormap('default');
229
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
230
caxis(c_map.caxis);
231
q=colorbar('vert');
232
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);
233
 
234
% Add the grid and the coast lines to the plot
235
m_grid;
236
m_coast('linewidth',1,'color','k');
237
title(num2str(z_uu.aklay(ilev)));
238
 
239
% Save figure
240
pre='';
241
if ( z_uu.aklay(ilev) < 10 )
242
  pre='0';
243
end
244
if ( z_uu.aklay(ilev) < 100 )
245
  pre=[pre '0' ];
246
end
247
if ( z_uu.aklay(ilev) < 1000 )
248
  pre=[pre '0' ];
249
end
250
if ( z_uu.aklay(ilev) < 10000 )
251
  pre=[pre '0' ];
252
end
253
 
254
figname = [ base '/ua_2d_' filename '_' pre num2str(z_uu.aklay(ilev)) '.eps' ];
255
set(gcf, 'PaperPosition', [2 1 15 10]);
256
print('-depsc2','-r0',figname);
257
 
258
% End loop over all levels
259
end