Subversion Repositories pvinversion.ecmwf

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed


folder  ='/lhome/sprenger/PV_Inversion_Tool/';
filename='Z1_20060115_18';

% -------------------------------------------------------------------------- 
% Plot qgPV anomaly (horizontal sections)
% --------------------------------------------------------------------------

% Base directory and filename
base = '/net/rossby/lhome/sprenger/PV_Inversion_Tool/';
folder   = base;
filename = 'Z1_20060115_18';
disp([folder filename])

% First image (otherwise first image is not correctly written)
figname = [base '/test.eps'];
close;
fh=figure('Units','pixels','Position',[100 100 900 900])
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);

% Load variables from file (on model levels)
m_pv = cdf_loadV(folder,filename,'QGPV_ANOM');

% Loop over all levels
for ilev=1:m_pv.nz

% Create a new figure
close;
fh=figure('Units','pixels','Position',[100 100 900 900])

% Set the geographical projection
m_proj('Equidistant Cylindrical','long',[m_pv.lonmin m_pv.lonmax],'lat',[m_pv.latmin m_pv.latmax]);

% Scale the plotting field for color map
fld=1e6*squeeze(m_pv.var(ilev,:,:));
c_map = scale_col(0:50:600,fld);

% Plot PV
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
[C,h]=m_contourf(lon,lat,c_map.data,c_map.xtick);
for icnt = 1: length(h)
  set( h(icnt), 'EdgeColor', 'none' )
end

% Add color bar
colormap('default');
ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,1);
caxis(c_map.caxis);
q=colorbar('hori');
set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label);

% Add the grid and the coast lines to the plot
m_grid;
m_coast('linewidth',1,'color','k');
title(num2str(m_pv.aklay(ilev)));

% Save figure
pre='';
if ( m_pv.aklay(ilev) < 10 )
  pre='0';
end
if ( m_pv.aklay(ilev) < 100 )
  pre=[pre '0' ];
end
if ( m_pv.aklay(ilev) < 1000 )
  pre=[pre '0' ];
end
if ( m_pv.aklay(ilev) < 10000 )
  pre=[pre '0' ];
end

figname = [ base '/qg_2d_' filename '_' pre num2str(m_pv.aklay(ilev)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);

% End loop over all levels
end