Subversion Repositories pvinversion.ecmwf

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

% 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)
z_th = cdf_loadV(folder,filename,'TH');
z_uu = cdf_loadV(folder,filename,'U');
z_vv = cdf_loadV(folder,filename,'V');

% -------------------------------------------------------------------------
% Lower boundary condition for potential temperature
% -------------------------------------------------------------------------

ilev=1;

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

% Set projection
load coast
h=axesm('MapProjection','eqdcylin');
setm(gca,'FLatLimit',[z_th.latmin z_th.latmax],'FLonLimit',[z_th.lonmin z_th.lonmax]);
h=plotm(lat,long);
gridm;

% Scale the plotting field for color map
fld=squeeze(z_th.var(ilev,:,:));
c_map = scale_col(230:5:310,fld);

% Plot TH
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]=contourfm(lat,lon,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,2);
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
title([ 'Lower Boundary Condition : Potential Temperature @' num2str(z_th.aklay(ilev)) ' m ASL' ]);

% Save figure
figname = [ base '/bound_lower_th_' filename '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);

% -------------------------------------------------------------------------
% Upper boundary condition for potential temperature
% -------------------------------------------------------------------------

ilev=z_th.nz;

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

% Set projection
load coast
h=axesm('MapProjection','eqdcylin');
setm(gca,'FLatLimit',[z_th.latmin z_th.latmax],'FLonLimit',[z_th.lonmin z_th.lonmax]);
h=plotm(lat,long);
gridm;

% Scale the plotting field for color map
fld=squeeze(z_th.var(ilev,:,:));
c_map = scale_col(450:10:550,fld);

% Plot TH
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]=contourfm(lat,lon,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,2);
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
title([ 'Upper Boundary Condition : Potential Temperature @ ' num2str(z_th.aklay(ilev)) ' m ASL' ]);

% Save figure
figname = [ base '/bound_upper_th_' filename '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);


% -------------------------------------------------------------------------
% Southern lateral boundary condition for zonal wind
% -------------------------------------------------------------------------

ilev=1;

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

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

% Plot
lev=z_uu.aklev/10000;
lon=z_uu.xmin + (0:z_uu.nx-1) * z_uu.dx;
[C,h]=contourf(lon,lev,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,2);
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
title([ 'Southern lateral Boundary Condition : Zonal Wind @ ' num2str(z_uu.ymin) '^\circ N' ]);
xlabel('Longitude');
ylabel('Height [km]');

% Save figure
figname = [ base '/bound_south_uu_' filename '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);


% -------------------------------------------------------------------------
% Northern lateral boundary condition for zonal wind
% -------------------------------------------------------------------------

ilev=z_uu.ny;

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

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

% Plot
lev=z_uu.aklev/10000;
lon=z_uu.xmin + (0:z_uu.nx-1) * z_uu.dx;
[C,h]=contourf(lon,lev,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,2);
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
title([ 'Northern lateral Boundary Condition : Zonal Wind @ ' num2str(z_uu.ymax) '^\circ N' ]);
xlabel('Longitude');
ylabel('Height [km]');

% Save figure
figname = [ base '/bound_north_uu_' filename '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);

% -------------------------------------------------------------------------
% Western lateral boundary condition for meridional wind
% -------------------------------------------------------------------------

ilev=1;

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

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

% Plot
lev=z_vv.aklev/10000;
lat=z_vv.ymin + (0:z_vv.ny-1) * z_vv.dy;
[C,h]=contourf(lat,lev,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,2);
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
title([ 'Western lateral Boundary Condition : Meridional Wind @ ' num2str(z_vv.xmin) '^\circ E' ]);
xlabel('Latitude');
ylabel('Height [km]');

% Save figure
figname = [ base '/bound_west_vv_' filename '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);

% -------------------------------------------------------------------------
% Eastern lateral boundary condition for meridional wind
% -------------------------------------------------------------------------

ilev=z_vv.nx;

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

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

% Plot
lev=z_vv.aklev/10000;
lat=z_vv.ymin + (0:z_vv.ny-1) * z_vv.dy;
[C,h]=contourf(lat,lev,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,2);
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
title([ 'Eastern lateral Boundary Condition : Meridional Wind @ ' num2str(z_vv.xmax) '^\circ E' ]);
xlabel('Latitude');
ylabel('Height [km]');

% Save figure
figname = [ base '/bound_east_vv_' filename '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);