function [ data, i4, r4, nam, xgrid, ygrid, history ] = nc_read( filename, ...
  record, levelreq );

%function [ data, i4, r4, nam, xgrid, ygrid, history ] = nc_read( filename, ...
%  record, levelreq );
%
% Reads one map of netCDF-formatted data.
%
% You need to include the following library to use this script.
% path( path, '/home/disk/tao/willa/matlab-epic' );
% 
%
% data		a row vector of data (of length nlat*nlon). -1 for end-of-file.
% i4, r4	metadata vectors (see Alexis' format specifications).
% 		The only non-zero elements in these vectors are the
%               following:
% iyear		year [i4(1)]
% ipd           month or season [i4(2)]
% iday          day [i4(3)]
% ihour		hour [i4(4)]
% nlon, nlat    number of longitudes and latitudes, respectively.  [i4(6:7)]
% miss   	The missing value flag. [i4(14)]  ***** For now this is
%               always NaN. *****
% base,fctr	add_offset and scale_factor for packing the data,
%               respectively. [r4(1:2)]
%               The use of these variables to pack and unpack the data in
%               the "CDC netCDF Standard Format" specifications and Alexis' 
%               pkd format are consistent.
% xgrid, ygrid  vectors of the longitude and latitude gridpoints,
%               respectively.
% history       an extended header (real words describing the data set
%               or calculation), if it exists
%
% filename      filename of the file to be read.
% record        number of the record to be read (default value 1)
% levelreq      (optional) vertical level to read data for, e.g., 1000
%               for 1000mb.  Type "ncdump -v level filename.nc" to see
%               how the levels are named.  The input is a number so no quotes.
%
% At present the script has no way to read a file in which multiple
% variables are given different variable names, for example, 
% geopotential height and temperature.
%
% July 1995: original script
%
% July 1996: I put in code to read in data which are intrinsically
% two-dimensional, for example, elevations.  If there are only two dimensions,
% the algorithm will fill output time information as zeros.  This is an
% ad hoc solution.  Obviously, you would also want this routine to read
% in time-latitude and time-longitude sections correctly.  
%
% November 1996: The time variable can be written one of two ways:
% 1) "so many units of time since some reference time," (COARDS standard),
% 2) "yyyymmddhhmmss", which is the older convention.
%
% 26 November 1996: Can now pick off the data for a selected level/height
% in a data set with multiple levels/heights.
%
% August 2000: Sets "scale_factor" and "add_offset" to 1 and 0, respectively
% if none are supplied.
%
% Todd Mitchell
%

if nargin==1 ; record = 1; end

i4 = zeros( 1, 30 );
r4 = zeros( 1, 20 );

%OPEN FILE
cdfid = mexcdf('OPEN',filename,'NOWRITE');
if cdfid<0; error('Cannot open file'); end

% GET AND CHECK NUMBER OF DIMENSIONS

[ ndim, nvar, natt, recdim, status ] = mexcdf( 'INQUIRE', cdfid );

while 0 
% Turn off the following check.  I have begun to write the time 
% variable twice: as hours since some reference time and 
% as yyyymmddhhmm, so the number of variables no longer bears the 
% simple relationship below.
if ndim~=nvar-1 & record==1
disp( sprintf( '\n ndim, nvar: %d %d', ndim, nvar ) );
disp( sprintf( 'Caution: ndim~=nvar-1 as expected.' ) );
end
end

if ndim<3 ; 
disp( sprintf( ...
  'Trouble: the number of dimensions, %d, is smaller than expected.', ndim ) );
error 
end

% Read in the attributes.
lennam  = zeros(nvar,1 );
natt    = zeros(nvar,1 );
ndim    = zeros(nvar,1 );
varnam  = zeros(nvar,20);
dimnam  = zeros(nvar,20);

for ivar = 1: nvar
[ nam, dtype, ndim(ivar), dimid, natt(ivar), status ] = ...
  mexcdf( 'VARINQ', cdfid, ivar-1 );
if ndim(ivar)>1; 
% This is the variable with the data of interest.
dimid_var=dimid; 
varid    =ivar; 
end
varnam(ivar,1:length(nam)) = nam;
end

% Determine the order in which the variables are written.
leng    = zeros(nvar,1);
for icnt = 1: length(dimid_var)
[ namdd, lend, status ] = mexcdf( 'DIMINQ', cdfid, dimid_var(icnt) );
for icnt2 = 1: nvar
if varnam(icnt2,1:length(namdd))==namdd; leng(icnt2) = lend; end
end
end  

% Read in all of the variables, but the 'data.'
% Test to see if the variable is the data.  Reading the data requires 
% knowledge of the size of the array, and so you have to read in the
% data last.  The if i~=varid & leng(i)>0 test keeps the algorithm
% from reading the data.
for ivar = 1: nvar
if ivar~=varid & leng(ivar)>0
% Read in time, depth, lat, or lon.
  d = mexcdf( 'VARGET', cdfid, ivar-1, 0, leng(ivar) ); 
  if varnam(ivar,1)=='t' | varnam(ivar,1)=='T'
    time = d;
    time_units_string = mexcdf( 'ATTGET', cdfid, ivar-1, 'units' );
% Test to see if the first letter of the variable name begins with
% an upper or lower case "x", or an upper or lower case "o" (as in
% longitude).  If yes, than this variable is the x-coordinate (the
% longitude).
  elseif varnam(ivar,1)=='x' | varnam(ivar,1)=='X' | varnam(ivar,2)=='o' ...
    | varnam(ivar,2)=='O'
    xgrid=d;
  elseif varnam(ivar,1)=='y' | varnam(ivar,1)=='Y' | varnam(ivar,2)=='a' ...
    | varnam(ivar,2)=='A'
    ygrid=d;
  else
    level = d;
  end
%disp( sprintf( 'Variable %d size %d %d', ivar, size(d) ) );
end
end

% Special case: Requested a map beyond the end of the data.  For example,
% a file of daily data for the present year where the number of maps
% is equal to the julian date.
if record>length(time) 
disp( sprintf( 'End-of-file.  "data" is set to -1.' ) )
data = -1;
break
end

% Read in the scale_factor and add_offset.
[ datatype, len, status ] = ...
  mexcdf( 'ATTINQ', cdfid, varid-1, 'scale_factor' );
if status==-1
  scale_factor = 1;
else
  scale_factor = mexcdf( 'ATTGET', cdfid, varid-1, 'scale_factor' );
end
[ datatype, len, status ] = ...
  mexcdf( 'ATTINQ', cdfid, varid-1, 'add_offset' );
if status==-1
  add_offset = 0;
else
  add_offset   = mexcdf( 'ATTGET', cdfid, varid-1, 'add_offset'   );
end
nam          = mexcdf( 'ATTGET', cdfid, varid-1, 'units'        );

% Read in the history variable
[ datatype, len, status ] = mexcdf( 'ATTINQ', cdfid, -1, 'history' );
if status~=-1
  history = mexcdf( 'ATTGET', cdfid, -1, 'history' );
else
  history = [];
end

% Read in the data variable.
% The first case is data that does not inherently have a time, such as
% topography.
if ndim<=3
[ data, status ] = mexcdf( 'VARGET', cdfid, varid-1, [record-1 0 0], ...
  [ 1 length(ygrid) length(xgrid) ] );
else
% Second case: single level data that changes with time.
  if nargin==2
  [ data, status ] = mexcdf( 'VARGET', cdfid, varid-1, ...
    [ record-1 0 0 0 ], [1 1 length(ygrid) length(xgrid) ] );
  else
% Third case: picking a single level off of a multi-level data set with
% time dependence.
  levelid = find( level==levelreq );
if isempty(levelid);
error( ...
  'You are not specifying the level correctly or this level does not exist.' );
end
disp( sprintf( 'levelid: %d', levelid ) )
  [ data, status ] = mexcdf( 'VARGET', cdfid, varid-1, ...
    [ record-1 levelid-1 0 0 ], [1 1 length(ygrid) length(xgrid) ] );
  end
end

% Convert the time variable into year, month, day, and hour.
% Two-dimensional data sets (for example, elevation) do not have a time
% value.
if ndim==2
iyear = 0;
ipd   = 1;
iday  = 1;
ihour = 0;
isec  = 0;
else

% There are at least two possibilities for the time convention:
% 1) time is written the old-fashioned way: yyyymmddhhmmss
% 2) time is written as "so many units since some reference time"

if isempty( findstr(time_units_string,'since') )
% Time is, by default, written as yyyymmddhhmmss
temp  = time(record);

else
% Time is written as "so many units since some reference time"
% Here, there are again two possibilities.  I wrote my reference times
% to include hours, minutes, and seconds so that it conforms with the
% files written at CDC.  Christian skips the hours, minutes, and seconds
% information if the basic unit of time is "hours".
if     strcmp(time_units_string(1:3),'sec') ; time_units = 'seconds' ;
elseif strcmp(time_units_string(1:3),'Sec') ; time_units = 'seconds'   ;
elseif strcmp(time_units_string(1:3),'hou') ; time_units = 'hours'   ;
elseif strcmp(time_units_string(1:3),'Hou') ; time_units = 'hours'   ;
elseif strcmp(time_units_string(1:3),'day') ; time_units = 'days'    ;
elseif strcmp(time_units_string(1:3),'Day') ; time_units = 'days'    ;
else   error( 'The time unit is not seconds, hours, or days.' );
end
ltu     = length(time_units_string);
yr1     = findstr( time_units_string, 'since' ) + 6;
dashpos = findstr( time_units_string, '-' );
yr2     = dashpos(1)-1;
mon1    = yr2 + 2;
mon2    = dashpos(2)-1;
dy1     = mon2 + 2;
if dy1==ltu | dy1+1==ltu
% Christian's code where "days" are the units of time, and there is
% no hours, minutes, and seconds in the reference time information.
dy2 = ltu;
basetime = ...
    sscanf( time_units_string(yr1:yr2),   '%i' ) * 1.0e10 ...
  + sscanf( time_units_string(mon1:mon2), '%i' ) * 1.0e8  ...
  + sscanf( time_units_string(dy1:dy2),   '%i' ) * 1.0e6;
else
% Reference time includes hours, minutes, seconds information.
temp2   = findstr( time_units_string(dy1+1:ltu), ' ' );
dy2     = dy1 + temp2(1) - 1;
hr1     = dy2 + 2;
colonpos = findstr( time_units_string, ':' );
hr2     = colonpos(1) - 1;
mn1     = hr2 + 2;
mn2     = colonpos(2) - 1;
se1     = mn2 + 2;
periodpos = findstr( time_units_string, '.' );
if isempty(periodpos) ; periodpos = ltu + 1; end
%if periodpos==[] ; periodpos = ltu + 1; end
basetime = ...
    sscanf( time_units_string(yr1:yr2),   '%i' ) * 1.0e10 ...
  + sscanf( time_units_string(mon1:mon2), '%i' ) * 1.0e8 ...
  + sscanf( time_units_string(dy1:dy2),   '%i' ) * 1.0e6 ...
  + sscanf( time_units_string(hr1:hr2),   '%i' ) * 1.0e4 ...
  + sscanf( time_units_string(mn1:mn2),   '%i' ) * 1.0e2 ...
  + sscanf( time_units_string(se1:periodpos-1), '%i' );
end

reltime  = time(record);
date = mexeps( 'array19tot', reltime, time_units, basetime );
temp = date;
end

iyear = fix( temp / 1.0e10 );
temp  = rem( temp,  1.0e10 );
ipd   = fix( temp / 1.0e8  );
temp  = rem( temp,  1.0e8  );
iday  = fix( temp / 1.0e6  );
temp  = rem( temp,  1.0e6  );
ihour = fix( temp / 1.0e4  );
temp  = rem( temp,  1.0e4  );
imin  = fix( temp / 1.0e2  );
temp  = rem( temp,  1.0e2  );
isec  = fix( temp ) + imin * 60;
end

% Put in NaNs for the missing values.
missing_value = mexcdf( 'ATTGET', cdfid, varid-1, 'missing_value' );
NaN2 = NaN;
% The CDC folks use INF as their missing value.  You need to 
% test for this slightly differently than you would other missing
% values.
if missing_value==Inf
x = find( data==Inf );
elseif missing_value==32767;
x = find( data==32767 );
else
x = find( abs(data(:)-missing_value) < 1.0e-4 );
end

if( length(x)>0 ); data(x) = ones( length(x), 1 ) * NaN2; end
NaN  = NaN2;
temp = find( data==32767 );

data = ( data * scale_factor ) + add_offset;

% Put the metadata into the "i4" and "r4" vectors.
nlon     = length(xgrid);
nlat     = length(ygrid);
base     = add_offset;
fctr     = scale_factor;
x0       = xgrid(1);
y0       = ygrid(1);
dx       = xgrid(2) - xgrid(1);
dy       = ygrid(2) - ygrid(1);
i4(1:5)  = [ iyear ipd iday ihour isec ];
i4(6:7)  = [ nlon nlat ];
i4(14)   = NaN;
r4(1:2)  = [ base fctr ];
r4(9:12) = [ x0 y0 dx dy ];

% Reshape the map of data so that it is a row vector (so that the 
% data is a more typical MATLAB form).
% The first transform puts the data in "map-format" --- it is arranged
% like a map.  The second transform operation puts the data into the a 
% row in the way that MATLAB expects things to be ordered.
data = data' ;
data = reshape( data', 1, nlon*nlat );

varid=varid-1;

mexcdf( 'CLOSE', cdfid );

disp( sprintf( ...
  'Finished reading record %d .  Year, month/season, day %d %d %d', ...
  record, iyear, ipd, iday ) )