Subversion Repositories pvinversion.ecmwf

% -------------------------------------------------------------------------
% Load files
% -------------------------------------------------------------------------

% Base directory and filename
base = '/net/rossby/lhome/sprenger/PV_Inversion_Tool/run/';
folder   = base;
filename = 'ANO_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)
z_th   = cdf_loadV(folder,filename,'TH');
z_uu   = cdf_loadV(folder,filename,'U');
z_vv   = cdf_loadV(folder,filename,'V');
z_qgpv = cdf_loadV(folder,filename,'QGPV');
z_psi  = cdf_loadV(folder,filename,'PSI');

% -------------------------------------------------------------------------
% Plot stream function
% -------------------------------------------------------------------------

for ilev=1:z_psi.nz

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

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

% Scale the plotting field for color map
fld=squeeze(z_psi.var(ilev,:,:));
c_map = scale_col(-4000:500:4000,fld);

% Plot Stream Function PSI
lat=z_psi.ymin + (0:z_psi.ny-1) * z_psi.dy;
lon=z_psi.xmin + (0:z_psi.nx-1) * z_psi.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('vert');
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);

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

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

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

% End loop over all levels
end

% -------------------------------------------------------------------------
% Plot anomaly of potential temperature
% -------------------------------------------------------------------------

for ilev=1:z_th.nz

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

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

% Scale the plotting field for color map
fld=squeeze(z_th.var(ilev,:,:));
c_map = scale_col(-20:2.5:20,fld);

% Plot potential temperature
lat=z_th.ymin + (0:z_th.ny-1) * z_th.dy;
lon=z_th.xmin + (0:z_th.nx-1) * z_th.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('vert');
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);

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

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

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

% End loop over all levels
end

% -------------------------------------------------------------------------
% Plot anomaly of meridional wind
% -------------------------------------------------------------------------

for ilev=1:z_vv.nz

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

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

% Scale the plotting field for color map
fld=squeeze(z_vv.var(ilev,:,:));
c_map = scale_col(-30:3:30,fld);

% Plot potential temperature
lat=z_vv.ymin + (0:z_vv.ny-1) * z_vv.dy;
lon=z_vv.xmin + (0:z_vv.nx-1) * z_vv.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('vert');
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);

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

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

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

% End loop over all levels
end

% -------------------------------------------------------------------------
% Plot anomaly of zonal wind
% -------------------------------------------------------------------------

for ilev=1:z_uu.nz

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

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

% Scale the plotting field for color map
fld=squeeze(z_uu.var(ilev,:,:));
c_map = scale_col(-30:3:30,fld);

% Plot potential temperature
lat=z_uu.ymin + (0:z_uu.ny-1) * z_uu.dy;
lon=z_uu.xmin + (0:z_uu.nx-1) * z_uu.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('vert');
set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label);

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

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

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

% End loop over all levels
end