Subversion Repositories pvinversion.ecmwf

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
 
2
% -------------------------------------------------------------------------
3
% Load files
4
% -------------------------------------------------------------------------
5
 
6
% Base directory and filename
7
base = '/net/rossby/lhome/sprenger/PV_Inversion_Tool/real/inp/';
8
folder   = base;
9
filename = 'H20060116_18';
10
disp([folder filename])
11
 
12
% First image (otherwise first image is not correctly written)
13
figname = [base '/test.eps'];
14
close;
15
fh=figure('Units','pixels','Position',[100 100 900 900])
16
set(gcf, 'PaperPosition', [2 1 15 10]);
17
print('-depsc2','-r0',figname);
18
 
19
% Load variables from file (on model levels)
20
m_th = cdf_loadV(folder,filename,'TH');
21
m_pv = cdf_loadV(folder,filename,'PV');
22
m_fl = cdf_loadV(folder,filename,'PV_FILT');
23
m_an = cdf_loadV(folder,filename,'PV_ANOM');
24
 
25
% Specify the plot domain
26
lonmin =-100;
27
lonmax = -20;
28
latmin =  10;
29
latmax =  80;
30
 
31
% Region for smoothing
32
lon1=-75;lon2=-50;
33
lat1= 30;lat2=55;
34
lev1=2000;lev2=10000;
35
 
36
% Plotting parameters
37
plotting=cellstr([ 'ew'; 'ns'; '2d'; 'pv'; 'an'; 'fl' ]);
38
%plotting=cellstr([ 'ew';  'pv'; 'an'; 'fl' ]);
39
 
40
% -------------------------------------------------------------------------
41
% PV (HORIZONTAL)
42
% -------------------------------------------------------------------------
43
 
44
if ( any(strcmp('2d',plotting)) & any(strcmp('pv',plotting)) )
45
 
46
disp('2D PV');
47
 
48
% Loop over all levels
49
for ilev=1:m_th.nz
50
 
51
% Create a new figure
52
close;
53
fh=figure('Units','pixels','Position',[100 100 900 900])
54
 
55
% Set the geographical projection
56
m_proj('Equidistant Cylindrical','long',[lonmin lonmax],'lat',[latmin latmax]);
57
 
58
% Scale the plotting field for color map
59
fld=squeeze(m_pv.var(ilev,:,:));
60
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
61
 
62
% Plot PV
63
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
64
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
65
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
66
for icnt = 1: length(h)
67
  set( h(icnt), 'EdgeColor', 'none' )
68
end
69
 
70
% Add color bar
71
colormap('default');
72
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
73
caxis(c_map.caxis);
74
q=colorbar('hori');
75
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
76
 
77
% Add the grid and the coast lines to the plot
78
m_grid;
79
m_coast('linewidth',1,'color','k');
80
title(num2str(m_th.aklay(ilev)));
81
 
82
% Add region of smoothing
83
if ( (m_th.aklay(ilev) > lev1) & (m_th.aklay(ilev) < lev2) ) 
84
  m_line([lon1 lon1],[lat1 lat2],'color','k','linewidth',2);
85
  m_line([lon1 lon2],[lat2 lat2],'color','k','linewidth',2);
86
  m_line([lon2 lon2],[lat2 lat1],'color','k','linewidth',2);
87
  m_line([lon1 lon2],[lat1 lat1],'color','k','linewidth',2);
88
end
89
 
90
% Save figure
91
pre='';
92
if ( m_th.aklay(ilev) < 10 )
93
  pre='0';
94
end
95
if ( m_th.aklay(ilev) < 100 )
96
  pre=[pre '0' ];
97
end
98
if ( m_th.aklay(ilev) < 1000 )
99
  pre=[pre '0' ];
100
end
101
if ( m_th.aklay(ilev) < 10000 )
102
  pre=[pre '0' ];
103
end
104
 
105
figname = [ base '/pv_2d_' filename '_' pre num2str(m_th.aklay(ilev)) '.eps' ];
106
set(gcf, 'PaperPosition', [2 1 15 10]);
107
print('-depsc2','-r0',figname);
108
 
109
% End loop over all levels
110
end
111
 
112
end
113
 
114
% -------------------------------------------------------------------------
115
% PV_FILT (HORIZONTAL)
116
% -------------------------------------------------------------------------
117
 
118
if ( any(strcmp('2d',plotting)) & any(strcmp('fl',plotting)) )
119
 
120
disp('2D FL');
121
 
122
% Loop over all levels
123
for ilev=1:m_th.nz
124
 
125
% Create a new figure
126
close;
127
fh=figure('Units','pixels','Position',[100 100 900 900])
128
 
129
% Set the geographical projection
130
m_proj('Equidistant Cylindrical','long',[lonmin lonmax],'lat',[latmin latmax]);
131
 
132
% Scale the plotting field for color map
133
fld=squeeze(m_fl.var(ilev,:,:));
134
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
135
 
136
% Plot PV
137
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
138
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
139
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
140
for icnt = 1: length(h)
141
  set( h(icnt), 'EdgeColor', 'none' )
142
end
143
 
144
% Add color bar
145
colormap('default');
146
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
147
caxis(c_map.caxis);
148
q=colorbar('hori');
149
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
150
 
151
% Add the grid and the coast lines to the plot
152
m_grid;
153
m_coast('linewidth',1,'color','k');
154
title(num2str(m_th.aklay(ilev)));
155
 
156
% Add region for smootrhing
157
if ( (m_th.aklay(ilev) > lev1) & (m_th.aklay(ilev) < lev2) ) 
158
  m_line([lon1 lon1],[lat1 lat2],'color','k','linewidth',2);
159
  m_line([lon1 lon2],[lat2 lat2],'color','k','linewidth',2);
160
  m_line([lon2 lon2],[lat2 lat1],'color','k','linewidth',2);
161
  m_line([lon1 lon2],[lat1 lat1],'color','k','linewidth',2);
162
end
163
 
164
% Save figure
165
pre='';
166
if ( m_th.aklay(ilev) < 10 )
167
  pre='0';
168
end
169
if ( m_th.aklay(ilev) < 100 )
170
  pre=[pre '0' ];
171
end
172
if ( m_th.aklay(ilev) < 1000 )
173
  pre=[pre '0' ];
174
end
175
if ( m_th.aklay(ilev) < 10000 )
176
  pre=[pre '0' ];
177
end
178
 
179
figname = [ base '/fl_2d_' filename '_' pre num2str(m_th.aklay(ilev)) '.eps' ];
180
set(gcf, 'PaperPosition', [2 1 15 10]);
181
print('-depsc2','-r0',figname);
182
 
183
% End loop over all levels
184
end
185
 
186
end
187
 
188
% -------------------------------------------------------------------------
189
% PV_ANOM (HORIZONTAL)
190
% -------------------------------------------------------------------------
191
 
192
if ( any(strcmp('2d',plotting)) & any(strcmp('an',plotting)) )
193
 
194
disp('2D AN');
195
 
196
% Loop over all levels
197
for ilev=1:m_th.nz
198
 
199
% Create a new figure
200
close;
201
fh=figure('Units','pixels','Position',[100 100 900 900])
202
 
203
% Set the geographical projection
204
m_proj('Equidistant Cylindrical','long',[lonmin lonmax],'lat',[latmin latmax]);
205
 
206
% Scale the plotting field for color map
207
fld=squeeze(m_an.var(ilev,:,:));
208
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
209
 
210
% Plot PV
211
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
212
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
213
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
214
for icnt = 1: length(h)
215
  set( h(icnt), 'EdgeColor', 'none' )
216
end
217
 
218
% Add color bar
219
colormap('default');
220
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
221
caxis(c_map.caxis);
222
q=colorbar('hori');
223
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
224
 
225
% Add the grid and the coast lines to the plot
226
m_grid;
227
m_coast('linewidth',1,'color','k');
228
title(num2str(m_th.aklay(ilev)));
229
 
230
% Add region for smootrhing
231
if ( (m_th.aklay(ilev) > lev1) & (m_th.aklay(ilev) < lev2) ) 
232
  m_line([lon1 lon1],[lat1 lat2],'color','k','linewidth',2);
233
  m_line([lon1 lon2],[lat2 lat2],'color','k','linewidth',2);
234
  m_line([lon2 lon2],[lat2 lat1],'color','k','linewidth',2);
235
  m_line([lon1 lon2],[lat1 lat1],'color','k','linewidth',2);
236
end
237
 
238
% Save figure
239
pre='';
240
if ( m_th.aklay(ilev) < 10 )
241
  pre='0';
242
end
243
if ( m_th.aklay(ilev) < 100 )
244
  pre=[pre '0' ];
245
end
246
if ( m_th.aklay(ilev) < 1000 )
247
  pre=[pre '0' ];
248
end
249
if ( m_th.aklay(ilev) < 10000 )
250
  pre=[pre '0' ];
251
end
252
 
253
figname = [ base '/an_2d_' filename '_' pre num2str(m_th.aklay(ilev)) '.eps' ];
254
set(gcf, 'PaperPosition', [2 1 15 10]);
255
print('-depsc2','-r0',figname);
256
 
257
% End loop over all levels
258
end
259
 
260
end
261
 
262
% -------------------------------------------------------------------------
263
% PV (VERTICAL EAST/WEST)
264
% -------------------------------------------------------------------------
265
 
266
if ( any(strcmp('ew',plotting)) & any(strcmp('pv',plotting)) )
267
 
268
disp('EW PV');
269
 
270
% Loop over all levels
271
for lat=latmin:latmax
272
 
273
% Get the y index of the latitude in the grid  
274
ilat=floor((lat-m_pv.ymin)/(m_pv.dy)+1); 
275
 
276
% Create a new figure
277
close;
278
fh=figure('Units','pixels','Position',[100 100 900 900])
279
 
280
% Scale the plotting field for color map
281
fld=squeeze(m_pv.var(:,ilat,:));
282
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
283
 
284
% Plot PV
285
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
286
lev=m_pv.aklay;
287
mask=(lon>=lonmin) & (lon<=lonmax);
288
[C,h]=contourf(lon(mask),lev,c_map.data(:,mask),c_map.xtick);
289
for icnt = 1: length(h)
290
  set( h(icnt), 'EdgeColor', 'none' )
291
end
292
 
293
% Add color bar
294
colormap('default');
295
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
296
caxis(c_map.caxis);
297
q=colorbar('hori');
298
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
299
 
300
% Add the grid and the coast lines to the plot
301
title([ 'Latitude = ' num2str(lat) ]);
302
 
303
% Add region for smootrhing
304
if ( (lat > lat1) & (lat < lat2) ) 
305
  line([lon1 lon2],[lev1 lev1],'color','k','linewidth',2);
306
  line([lon2 lon2],[lev1 lev2],'color','k','linewidth',2);
307
  line([lon2 lon1],[lev2 lev2],'color','k','linewidth',2);
308
  line([lon1 lon1],[lev2 lev1],'color','k','linewidth',2);
309
end
310
 
311
% Save figure
312
if ( lat < 0 ) 
313
  pre='s'
314
else
315
  pre='n'
316
end
317
if ( abs(lat) < 10 ) 
318
  pre=[ pre '0'];
319
end
320
 
321
figname = [ base '/pv_ew_' filename '_' pre num2str(abs(lat)) '.eps' ];
322
set(gcf, 'PaperPosition', [2 1 15 10]);
323
print('-depsc2','-r0',figname);
324
 
325
% End loop over all levels
326
end
327
 
328
end
329
 
330
% -------------------------------------------------------------------------
331
% PV_FILT (VERTICAL EAST/WEST)
332
% -------------------------------------------------------------------------
333
 
334
if ( any(strcmp('ew',plotting)) & any(strcmp('fl',plotting)) )
335
 
336
disp('EW FL');
337
 
338
% Loop over all levels
339
for lat=latmin:latmax
340
 
341
% Get the y index of the latitude in the grid  
342
ilat=floor((lat-m_pv.ymin)/(m_pv.dy)+1); 
343
 
344
% Create a new figure
345
close;
346
fh=figure('Units','pixels','Position',[100 100 900 900])
347
 
348
% Scale the plotting field for color map
349
fld=squeeze(m_fl.var(:,ilat,:));
350
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
351
 
352
% Plot PV_FILT
353
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
354
lev=m_pv.aklay;
355
mask=(lon>=lonmin) & (lon<=lonmax);
356
[C,h]=contourf(lon(mask),lev,c_map.data(:,mask),c_map.xtick);
357
for icnt = 1: length(h)
358
  set( h(icnt), 'EdgeColor', 'none' )
359
end
360
 
361
% Add color bar
362
colormap('default');
363
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
364
caxis(c_map.caxis);
365
q=colorbar('hori');
366
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
367
 
368
% Add the grid and the coast lines to the plot
369
title([ 'Latitude = ' num2str(lat) ]);
370
 
371
% Add region for smootrhing
372
if ( (lat > lat1) & (lat < lat2) ) 
373
  line([lon1 lon2],[lev1 lev1],'color','k','linewidth',2);
374
  line([lon2 lon2],[lev1 lev2],'color','k','linewidth',2);
375
  line([lon2 lon1],[lev2 lev2],'color','k','linewidth',2);
376
  line([lon1 lon1],[lev2 lev1],'color','k','linewidth',2);
377
end
378
 
379
% Save figure
380
if ( lat < 0 ) 
381
  pre='s'
382
else
383
  pre='n'
384
end
385
if ( abs(lat) < 10 ) 
386
  pre=[ pre '0'];
387
end
388
 
389
figname = [ base '/fl_ew_' filename '_' pre num2str(abs(lat)) '.eps' ];
390
set(gcf, 'PaperPosition', [2 1 15 10]);
391
print('-depsc2','-r0',figname);
392
 
393
% End loop over all levels
394
end
395
 
396
end
397
 
398
% -------------------------------------------------------------------------
399
% PV_ANOM (VERTICAL EAST/WEST)
400
% -------------------------------------------------------------------------
401
 
402
if ( any(strcmp('ew',plotting)) & any(strcmp('an',plotting)) )
403
 
404
disp('EW AN');
405
 
406
% Loop over all levels
407
for lat=latmin:latmax
408
 
409
% Get the y index of the latitude in the grid  
410
ilat=floor((lat-m_pv.ymin)/(m_pv.dy)+1); 
411
 
412
% Create a new figure
413
close;
414
fh=figure('Units','pixels','Position',[100 100 900 900])
415
 
416
% Scale the plotting field for color map
417
fld=squeeze(m_an.var(:,ilat,:));
418
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
419
 
420
% Plot PV_FILT
421
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
422
lev=m_pv.aklay;
423
mask=(lon>=lonmin) & (lon<=lonmax);
424
[C,h]=contourf(lon(mask),lev,c_map.data(:,mask),c_map.xtick);
425
for icnt = 1: length(h)
426
  set( h(icnt), 'EdgeColor', 'none' )
427
end
428
 
429
% Add color bar
430
colormap('default');
431
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
432
caxis(c_map.caxis);
433
q=colorbar('hori');
434
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
435
 
436
% Add the grid and the coast lines to the plot
437
title([ 'Latitude = ' num2str(lat) ]);
438
 
439
% Add region for smootrhing
440
if ( (lat > lat1) & (lat < lat2) ) 
441
  line([lon1 lon2],[lev1 lev1],'color','k','linewidth',2);
442
  line([lon2 lon2],[lev1 lev2],'color','k','linewidth',2);
443
  line([lon2 lon1],[lev2 lev2],'color','k','linewidth',2);
444
  line([lon1 lon1],[lev2 lev1],'color','k','linewidth',2);
445
end
446
 
447
% Save figure
448
if ( lat < 0 ) 
449
  pre='s'
450
else
451
  pre='n'
452
end
453
if ( abs(lat) < 10 ) 
454
  pre=[ pre '0'];
455
end
456
 
457
figname = [ base '/an_ew_' filename '_' pre num2str(abs(lat)) '.eps' ];
458
set(gcf, 'PaperPosition', [2 1 15 10]);
459
print('-depsc2','-r0',figname);
460
 
461
% End loop over all levels
462
end
463
 
464
end
465
 
466
% -------------------------------------------------------------------------
467
% PV (VERTICAL NORTH/SOUTH)
468
% -------------------------------------------------------------------------
469
 
470
if ( any(strcmp('ns',plotting)) & any(strcmp('pv',plotting)) )
471
 
472
disp('NS PV');
473
 
474
% Loop over all levels
475
for lon=lonmin:lonmax
476
 
477
% Get the y index of the latitude in the grid  
478
ilon=floor((lon-m_pv.xmin)/(m_pv.dx)+1); 
479
 
480
% Create a new figure
481
close;
482
fh=figure('Units','pixels','Position',[100 100 900 900])
483
 
484
% Scale the plotting field for color map
485
fld=squeeze(m_pv.var(:,:,ilon));
486
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
487
 
488
% Plot PV
489
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
490
lev=m_pv.aklay;
491
mask=(lat>=latmin) & (lat<=latmax);
492
[C,h]=contourf(lat(mask),lev,c_map.data(:,mask),c_map.xtick);
493
for icnt = 1: length(h)
494
  set( h(icnt), 'EdgeColor', 'none' )
495
end
496
 
497
% Add color bar
498
colormap('default');
499
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
500
caxis(c_map.caxis);
501
q=colorbar('hori');
502
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
503
 
504
% Add the grid and the coast lines to the plot
505
title([ 'Longitude = ' num2str(lon) ]);
506
 
507
% Add region for smootrhing
508
if ( (lon > lon1) & (lon < lon2) ) 
509
  line([lat1 lat2],[lev1 lev1],'color','k','linewidth',2);
510
  line([lat2 lat2],[lev1 lev2],'color','k','linewidth',2);
511
  line([lat2 lat1],[lev2 lev2],'color','k','linewidth',2);
512
  line([lat1 lat1],[lev2 lev1],'color','k','linewidth',2);
513
end
514
 
515
% Save figure
516
if ( lon < 0 ) 
517
  pre='w'
518
else
519
  pre='e'
520
end
521
if ( abs(lon) < 10 ) 
522
  pre=[ pre '00'];
523
elseif ( abs(lon) < 100 )
524
  pre=[ pre '0'];
525
end
526
 
527
figname = [ base '/pv_ns_' filename '_' pre num2str(abs(lon)) '.eps' ];
528
set(gcf, 'PaperPosition', [2 1 15 10]);
529
print('-depsc2','-r0',figname);
530
 
531
% End loop over all levels
532
end
533
 
534
end
535
 
536
% -------------------------------------------------------------------------
537
% PV_FILT (VERTICAL NORTH/SOUTH)
538
% -------------------------------------------------------------------------
539
 
540
if ( any(strcmp('ns',plotting)) & any(strcmp('fl',plotting)) )
541
 
542
disp('NS FL');
543
 
544
% Loop over all levels
545
for lon=lonmin:lonmax
546
 
547
% Get the y index of the latitude in the grid  
548
ilon=floor((lon-m_pv.xmin)/(m_pv.dx)+1); 
549
 
550
% Create a new figure
551
close;
552
fh=figure('Units','pixels','Position',[100 100 900 900])
553
 
554
% Scale the plotting field for color map
555
fld=squeeze(m_fl.var(:,:,ilon));
556
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
557
 
558
% Plot PV
559
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
560
lev=m_pv.aklay;
561
mask=(lat>=latmin) & (lat<=latmax);
562
[C,h]=contourf(lat(mask),lev,c_map.data(:,mask),c_map.xtick);
563
for icnt = 1: length(h)
564
  set( h(icnt), 'EdgeColor', 'none' )
565
end
566
 
567
% Add color bar
568
colormap('default');
569
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
570
caxis(c_map.caxis);
571
q=colorbar('hori');
572
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
573
 
574
% Add the grid and the coast lines to the plot
575
title([ 'Longitude = ' num2str(lon) ]);
576
 
577
% Add region for smootrhing
578
if ( (lon > lon1) & (lon < lon2) ) 
579
  line([lat1 lat2],[lev1 lev1],'color','k','linewidth',2);
580
  line([lat2 lat2],[lev1 lev2],'color','k','linewidth',2);
581
  line([lat2 lat1],[lev2 lev2],'color','k','linewidth',2);
582
  line([lat1 lat1],[lev2 lev1],'color','k','linewidth',2);
583
end
584
 
585
% Save figure
586
if ( lon < 0 ) 
587
  pre='w'
588
else
589
  pre='e'
590
end
591
if ( abs(lon) < 10 ) 
592
  pre=[ pre '00'];
593
elseif ( abs(lon) < 100 )
594
  pre=[ pre '0'];
595
end
596
 
597
figname = [ base '/fl_ns_' filename '_' pre num2str(abs(lon)) '.eps' ];
598
set(gcf, 'PaperPosition', [2 1 15 10]);
599
print('-depsc2','-r0',figname);
600
 
601
% End loop over all levels
602
end
603
 
604
end
605
 
606
% -------------------------------------------------------------------------
607
% PV_ANOM (VERTICAL NORTH/SOUTH)
608
% -------------------------------------------------------------------------
609
 
610
if ( any(strcmp('ns',plotting)) & any(strcmp('an',plotting)) )
611
 
612
disp('NS AN');
613
 
614
% Loop over all levels
615
for lon=lonmin:lonmax
616
 
617
% Get the y index of the latitude in the grid  
618
ilon=floor((lon-m_pv.xmin)/(m_pv.dx)+1); 
619
 
620
% Create a new figure
621
close;
622
fh=figure('Units','pixels','Position',[100 100 900 900])
623
 
624
% Scale the plotting field for color map
625
fld=squeeze(m_an.var(:,:,ilon));
626
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
627
 
628
% Plot PV
629
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
630
lev=m_pv.aklay;
631
mask=(lat>=latmin) & (lat<=latmax);
632
[C,h]=contourf(lat(mask),lev,c_map.data(:,mask),c_map.xtick);
633
for icnt = 1: length(h)
634
  set( h(icnt), 'EdgeColor', 'none' )
635
end
636
 
637
% Add color bar
638
colormap('default');
639
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
640
caxis(c_map.caxis);
641
q=colorbar('hori');
642
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);
643
 
644
% Add the grid and the coast lines to the plot
645
title([ 'Longitude = ' num2str(lon) ]);
646
 
647
% Add region for smootrhing
648
if ( (lon > lon1) & (lon < lon2) ) 
649
  line([lat1 lat2],[lev1 lev1],'color','k','linewidth',2);
650
  line([lat2 lat2],[lev1 lev2],'color','k','linewidth',2);
651
  line([lat2 lat1],[lev2 lev2],'color','k','linewidth',2);
652
  line([lat1 lat1],[lev2 lev1],'color','k','linewidth',2);
653
end
654
 
655
% Save figure
656
if ( lon < 0 ) 
657
  pre='w'
658
else
659
  pre='e'
660
end
661
if ( abs(lon) < 10 ) 
662
  pre=[ pre '00'];
663
elseif ( abs(lon) < 100 )
664
  pre=[ pre '0'];
665
end
666
 
667
figname = [ base '/an_ns_' filename '_' pre num2str(abs(lon)) '.eps' ];
668
set(gcf, 'PaperPosition', [2 1 15 10]);
669
print('-depsc2','-r0',figname);
670
 
671
% End loop over all levels
672
end
673
 
674
end