Subversion Repositories pvinversion.ecmwf

Compare Revisions

No changes between revisions

Ignore whitespace Rev 2 → Rev 3

/trunk/prep/def_anomaly.m
0,0 → 1,674
 
% -------------------------------------------------------------------------
% Load files
% -------------------------------------------------------------------------
 
% Base directory and filename
base = '/net/rossby/lhome/sprenger/PV_Inversion_Tool/real/inp/';
folder = base;
filename = 'H20060116_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_th = cdf_loadV(folder,filename,'TH');
m_pv = cdf_loadV(folder,filename,'PV');
m_fl = cdf_loadV(folder,filename,'PV_FILT');
m_an = cdf_loadV(folder,filename,'PV_ANOM');
 
% Specify the plot domain
lonmin =-100;
lonmax = -20;
latmin = 10;
latmax = 80;
 
% Region for smoothing
lon1=-75;lon2=-50;
lat1= 30;lat2=55;
lev1=2000;lev2=10000;
 
% Plotting parameters
plotting=cellstr([ 'ew'; 'ns'; '2d'; 'pv'; 'an'; 'fl' ]);
%plotting=cellstr([ 'ew'; 'pv'; 'an'; 'fl' ]);
 
% -------------------------------------------------------------------------
% PV (HORIZONTAL)
% -------------------------------------------------------------------------
 
if ( any(strcmp('2d',plotting)) & any(strcmp('pv',plotting)) )
 
disp('2D PV');
% Loop over all levels
for ilev=1:m_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',[lonmin lonmax],'lat',[latmin latmax]);
 
% Scale the plotting field for color map
fld=squeeze(m_pv.var(ilev,:,:));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],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_th.aklay(ilev)));
 
% Add region of smoothing
if ( (m_th.aklay(ilev) > lev1) & (m_th.aklay(ilev) < lev2) )
m_line([lon1 lon1],[lat1 lat2],'color','k','linewidth',2);
m_line([lon1 lon2],[lat2 lat2],'color','k','linewidth',2);
m_line([lon2 lon2],[lat2 lat1],'color','k','linewidth',2);
m_line([lon1 lon2],[lat1 lat1],'color','k','linewidth',2);
end
 
% Save figure
pre='';
if ( m_th.aklay(ilev) < 10 )
pre='0';
end
if ( m_th.aklay(ilev) < 100 )
pre=[pre '0' ];
end
if ( m_th.aklay(ilev) < 1000 )
pre=[pre '0' ];
end
if ( m_th.aklay(ilev) < 10000 )
pre=[pre '0' ];
end
 
figname = [ base '/pv_2d_' filename '_' pre num2str(m_th.aklay(ilev)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV_FILT (HORIZONTAL)
% -------------------------------------------------------------------------
 
if ( any(strcmp('2d',plotting)) & any(strcmp('fl',plotting)) )
 
disp('2D FL');
 
% Loop over all levels
for ilev=1:m_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',[lonmin lonmax],'lat',[latmin latmax]);
 
% Scale the plotting field for color map
fld=squeeze(m_fl.var(ilev,:,:));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],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_th.aklay(ilev)));
 
% Add region for smootrhing
if ( (m_th.aklay(ilev) > lev1) & (m_th.aklay(ilev) < lev2) )
m_line([lon1 lon1],[lat1 lat2],'color','k','linewidth',2);
m_line([lon1 lon2],[lat2 lat2],'color','k','linewidth',2);
m_line([lon2 lon2],[lat2 lat1],'color','k','linewidth',2);
m_line([lon1 lon2],[lat1 lat1],'color','k','linewidth',2);
end
 
% Save figure
pre='';
if ( m_th.aklay(ilev) < 10 )
pre='0';
end
if ( m_th.aklay(ilev) < 100 )
pre=[pre '0' ];
end
if ( m_th.aklay(ilev) < 1000 )
pre=[pre '0' ];
end
if ( m_th.aklay(ilev) < 10000 )
pre=[pre '0' ];
end
 
figname = [ base '/fl_2d_' filename '_' pre num2str(m_th.aklay(ilev)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV_ANOM (HORIZONTAL)
% -------------------------------------------------------------------------
 
if ( any(strcmp('2d',plotting)) & any(strcmp('an',plotting)) )
 
disp('2D AN');
 
% Loop over all levels
for ilev=1:m_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',[lonmin lonmax],'lat',[latmin latmax]);
 
% Scale the plotting field for color map
fld=squeeze(m_an.var(ilev,:,:));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],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_th.aklay(ilev)));
 
% Add region for smootrhing
if ( (m_th.aklay(ilev) > lev1) & (m_th.aklay(ilev) < lev2) )
m_line([lon1 lon1],[lat1 lat2],'color','k','linewidth',2);
m_line([lon1 lon2],[lat2 lat2],'color','k','linewidth',2);
m_line([lon2 lon2],[lat2 lat1],'color','k','linewidth',2);
m_line([lon1 lon2],[lat1 lat1],'color','k','linewidth',2);
end
 
% Save figure
pre='';
if ( m_th.aklay(ilev) < 10 )
pre='0';
end
if ( m_th.aklay(ilev) < 100 )
pre=[pre '0' ];
end
if ( m_th.aklay(ilev) < 1000 )
pre=[pre '0' ];
end
if ( m_th.aklay(ilev) < 10000 )
pre=[pre '0' ];
end
 
figname = [ base '/an_2d_' filename '_' pre num2str(m_th.aklay(ilev)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV (VERTICAL EAST/WEST)
% -------------------------------------------------------------------------
 
if ( any(strcmp('ew',plotting)) & any(strcmp('pv',plotting)) )
disp('EW PV');
 
% Loop over all levels
for lat=latmin:latmax
 
% Get the y index of the latitude in the grid
ilat=floor((lat-m_pv.ymin)/(m_pv.dy)+1);
 
% Create a new figure
close;
fh=figure('Units','pixels','Position',[100 100 900 900])
 
% Scale the plotting field for color map
fld=squeeze(m_pv.var(:,ilat,:));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
 
% Plot PV
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
lev=m_pv.aklay;
mask=(lon>=lonmin) & (lon<=lonmax);
[C,h]=contourf(lon(mask),lev,c_map.data(:,mask),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
title([ 'Latitude = ' num2str(lat) ]);
 
% Add region for smootrhing
if ( (lat > lat1) & (lat < lat2) )
line([lon1 lon2],[lev1 lev1],'color','k','linewidth',2);
line([lon2 lon2],[lev1 lev2],'color','k','linewidth',2);
line([lon2 lon1],[lev2 lev2],'color','k','linewidth',2);
line([lon1 lon1],[lev2 lev1],'color','k','linewidth',2);
end
 
% Save figure
if ( lat < 0 )
pre='s'
else
pre='n'
end
if ( abs(lat) < 10 )
pre=[ pre '0'];
end
 
figname = [ base '/pv_ew_' filename '_' pre num2str(abs(lat)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV_FILT (VERTICAL EAST/WEST)
% -------------------------------------------------------------------------
 
if ( any(strcmp('ew',plotting)) & any(strcmp('fl',plotting)) )
disp('EW FL');
 
% Loop over all levels
for lat=latmin:latmax
 
% Get the y index of the latitude in the grid
ilat=floor((lat-m_pv.ymin)/(m_pv.dy)+1);
 
% Create a new figure
close;
fh=figure('Units','pixels','Position',[100 100 900 900])
 
% Scale the plotting field for color map
fld=squeeze(m_fl.var(:,ilat,:));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
 
% Plot PV_FILT
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
lev=m_pv.aklay;
mask=(lon>=lonmin) & (lon<=lonmax);
[C,h]=contourf(lon(mask),lev,c_map.data(:,mask),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
title([ 'Latitude = ' num2str(lat) ]);
 
% Add region for smootrhing
if ( (lat > lat1) & (lat < lat2) )
line([lon1 lon2],[lev1 lev1],'color','k','linewidth',2);
line([lon2 lon2],[lev1 lev2],'color','k','linewidth',2);
line([lon2 lon1],[lev2 lev2],'color','k','linewidth',2);
line([lon1 lon1],[lev2 lev1],'color','k','linewidth',2);
end
 
% Save figure
if ( lat < 0 )
pre='s'
else
pre='n'
end
if ( abs(lat) < 10 )
pre=[ pre '0'];
end
 
figname = [ base '/fl_ew_' filename '_' pre num2str(abs(lat)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV_ANOM (VERTICAL EAST/WEST)
% -------------------------------------------------------------------------
 
if ( any(strcmp('ew',plotting)) & any(strcmp('an',plotting)) )
 
disp('EW AN');
 
% Loop over all levels
for lat=latmin:latmax
 
% Get the y index of the latitude in the grid
ilat=floor((lat-m_pv.ymin)/(m_pv.dy)+1);
 
% Create a new figure
close;
fh=figure('Units','pixels','Position',[100 100 900 900])
 
% Scale the plotting field for color map
fld=squeeze(m_an.var(:,ilat,:));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
 
% Plot PV_FILT
lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx;
lev=m_pv.aklay;
mask=(lon>=lonmin) & (lon<=lonmax);
[C,h]=contourf(lon(mask),lev,c_map.data(:,mask),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
title([ 'Latitude = ' num2str(lat) ]);
 
% Add region for smootrhing
if ( (lat > lat1) & (lat < lat2) )
line([lon1 lon2],[lev1 lev1],'color','k','linewidth',2);
line([lon2 lon2],[lev1 lev2],'color','k','linewidth',2);
line([lon2 lon1],[lev2 lev2],'color','k','linewidth',2);
line([lon1 lon1],[lev2 lev1],'color','k','linewidth',2);
end
 
% Save figure
if ( lat < 0 )
pre='s'
else
pre='n'
end
if ( abs(lat) < 10 )
pre=[ pre '0'];
end
 
figname = [ base '/an_ew_' filename '_' pre num2str(abs(lat)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV (VERTICAL NORTH/SOUTH)
% -------------------------------------------------------------------------
 
if ( any(strcmp('ns',plotting)) & any(strcmp('pv',plotting)) )
disp('NS PV');
% Loop over all levels
for lon=lonmin:lonmax
 
% Get the y index of the latitude in the grid
ilon=floor((lon-m_pv.xmin)/(m_pv.dx)+1);
 
% Create a new figure
close;
fh=figure('Units','pixels','Position',[100 100 900 900])
 
% Scale the plotting field for color map
fld=squeeze(m_pv.var(:,:,ilon));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
 
% Plot PV
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
lev=m_pv.aklay;
mask=(lat>=latmin) & (lat<=latmax);
[C,h]=contourf(lat(mask),lev,c_map.data(:,mask),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
title([ 'Longitude = ' num2str(lon) ]);
 
% Add region for smootrhing
if ( (lon > lon1) & (lon < lon2) )
line([lat1 lat2],[lev1 lev1],'color','k','linewidth',2);
line([lat2 lat2],[lev1 lev2],'color','k','linewidth',2);
line([lat2 lat1],[lev2 lev2],'color','k','linewidth',2);
line([lat1 lat1],[lev2 lev1],'color','k','linewidth',2);
end
 
% Save figure
if ( lon < 0 )
pre='w'
else
pre='e'
end
if ( abs(lon) < 10 )
pre=[ pre '00'];
elseif ( abs(lon) < 100 )
pre=[ pre '0'];
end
 
figname = [ base '/pv_ns_' filename '_' pre num2str(abs(lon)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV_FILT (VERTICAL NORTH/SOUTH)
% -------------------------------------------------------------------------
 
if ( any(strcmp('ns',plotting)) & any(strcmp('fl',plotting)) )
disp('NS FL');
 
% Loop over all levels
for lon=lonmin:lonmax
 
% Get the y index of the latitude in the grid
ilon=floor((lon-m_pv.xmin)/(m_pv.dx)+1);
 
% Create a new figure
close;
fh=figure('Units','pixels','Position',[100 100 900 900])
 
% Scale the plotting field for color map
fld=squeeze(m_fl.var(:,:,ilon));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
 
% Plot PV
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
lev=m_pv.aklay;
mask=(lat>=latmin) & (lat<=latmax);
[C,h]=contourf(lat(mask),lev,c_map.data(:,mask),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
title([ 'Longitude = ' num2str(lon) ]);
 
% Add region for smootrhing
if ( (lon > lon1) & (lon < lon2) )
line([lat1 lat2],[lev1 lev1],'color','k','linewidth',2);
line([lat2 lat2],[lev1 lev2],'color','k','linewidth',2);
line([lat2 lat1],[lev2 lev2],'color','k','linewidth',2);
line([lat1 lat1],[lev2 lev1],'color','k','linewidth',2);
end
 
% Save figure
if ( lon < 0 )
pre='w'
else
pre='e'
end
if ( abs(lon) < 10 )
pre=[ pre '00'];
elseif ( abs(lon) < 100 )
pre=[ pre '0'];
end
 
figname = [ base '/fl_ns_' filename '_' pre num2str(abs(lon)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
 
% -------------------------------------------------------------------------
% PV_ANOM (VERTICAL NORTH/SOUTH)
% -------------------------------------------------------------------------
 
if ( any(strcmp('ns',plotting)) & any(strcmp('an',plotting)) )
disp('NS AN');
 
% Loop over all levels
for lon=lonmin:lonmax
 
% Get the y index of the latitude in the grid
ilon=floor((lon-m_pv.xmin)/(m_pv.dx)+1);
 
% Create a new figure
close;
fh=figure('Units','pixels','Position',[100 100 900 900])
 
% Scale the plotting field for color map
fld=squeeze(m_an.var(:,:,ilon));
c_map = scale_col([-.5,-.2,0.,.2,.5,1.,1.5,2.,3.,4.,6.,8.,10.],fld);
 
% Plot PV
lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy;
lev=m_pv.aklay;
mask=(lat>=latmin) & (lat<=latmax);
[C,h]=contourf(lat(mask),lev,c_map.data(:,mask),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
title([ 'Longitude = ' num2str(lon) ]);
 
% Add region for smootrhing
if ( (lon > lon1) & (lon < lon2) )
line([lat1 lat2],[lev1 lev1],'color','k','linewidth',2);
line([lat2 lat2],[lev1 lev2],'color','k','linewidth',2);
line([lat2 lat1],[lev2 lev2],'color','k','linewidth',2);
line([lat1 lat1],[lev2 lev1],'color','k','linewidth',2);
end
 
% Save figure
if ( lon < 0 )
pre='w'
else
pre='e'
end
if ( abs(lon) < 10 )
pre=[ pre '00'];
elseif ( abs(lon) < 100 )
pre=[ pre '0'];
end
 
figname = [ base '/an_ns_' filename '_' pre num2str(abs(lon)) '.eps' ];
set(gcf, 'PaperPosition', [2 1 15 10]);
print('-depsc2','-r0',figname);
 
% End loop over all levels
end
 
end
Property changes:
Added: svn:executable