% Name/location of the netCDF file
nc_file = '~/job/test_data/ERA40.nc';

% Time level and vertical levels to plot
time_lev = 1;
z_lev = 1;

% Latitude and longitude range to plot over
lat_range = [50 80];
lon_range = [300 360];

% Map projection
projection = 'lambert';

% Fields to contour/shade. If more than one is specified then the l2 norm
% of these fields will be shaded (e.g. specifying u and v would shade speed).
% These should correspond to the variable names in the netCDF file.
field_contour = [{'U10'} {'V10'}];

% Fields to use to plot vectors. There must be either 0 or 2 fields here;
field_vector = [{'U10'} {'V10'}];

% Define the names of the latitude and the longitude dimensions in the 
% netCDF file.
lat_name = 'latitude';
lon_name = 'longitude';

%================= All options set above this point ====================%

% Read the latitude and longitude variable from file
lat = double(ncread(nc_file,lat_name));
lon = double(ncread(nc_file,lon_name));

% Save the length of lat & lon
num_lat = length(lat);
num_lon = length(lon);

% Loop over all of the variables in field_contour and calculate the magnitude of 
% these fields.
field = 0;
num_fields = length(field_contour);
for(i=1:length(field_contour))
  field = (field) + double((ncread(nc_file,char(field_contour(i)),[time_lev z_lev 1 1]-1,[1 1 num_lat num_lon]).^num_fields));
end
field = nthroot(field,num_fields);

% Trim down the size of the latitude and longitude arrays to be slightly larger than 
% the specified domain
lon_index = find(lon<max(lon_range)&lon>min(lon_range));
lon_index = [min(lon_index)-1 lon_index' max(lon_index)+1];
lon_index = lon_index(lon_index>0);
lon_index = lon_index(lon_index<=num_lon);

lat_index = find(lat<max(lat_range)&lat>min(lat_range));
lat_index = [min(lat_index-1) lat_index' max(lat_index)+1];
lat_index = lat_index(lat_index>0);
lat_index = lat_index(lat_index<=num_lat);

% Trim down field to the size of the specified domain
field = field(lat_index,lon_index);

% Set up the mapping projection; Shade specified field; grid and add coastlines
m_proj(projection,'lats',[min(lat_range) max(lat_range)],'lons',[min(lon_range) max(lon_range)])
m_pcolor(lon(lon_index),lat(lat_index),field),shading flat;
hold on;
m_coast('color','k','linewidth',2)
m_grid('box','fancy','linestyle','none');


% Read and add vectors to the plot, checking that only two of these are specified
num_fields = length(field_vector);
for(i=1:num_fields)
  eval(['field' num2str(i) '= double(ncread(nc_file,char(field_vector(i)),[time_lev z_lev 1 1]-1,[1 1 num_lat num_lon]));']);
end
if(num_fields==2)
  [lon lat] = meshgrid(lon,lat);
  m_quiver(lon(lat_index,lon_index),lat(lat_index,lon_index),field1(lat_index,lon_index),field2(lat_index,lon_index),'k');
elseif(num_fields~=0)
  error('field_vector must contain 0 or 2 field names');
end