Subversion Repositories pvinversion.ecmwf

% -------------------------------------------------------------------------
% 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