function [gregorian_time, serial_time, gregorian_base, serial_base, ...	  sizem, serial_time_jd, serial_base_jd] = ...    timenc(file, time_var, bl_corner, tr_corner, calendar)% TIMENC   Returns time information for a file that meets COARDS standards.%% function [gregorian_time, serial_time, gregorian_base, serial_base, ...% 	  sizem, serial_time_jd, serial_base_jd] = ...%         timenc(file, time_var, bl_corner, tr_corner, calendar);%% DESCRIPTION:%% timenc finds the time vector and the corresponding base date for a% netcdf file or DODS/OPEnDAP dataset that meets COARDS standards. In% practice this means that time-like variable should have a units attribute% of a certain form. An example is:%      'seconds since 1992-10-8 15:15:42.5 -6:00'% This indicates seconds since October 8th, 1992 at 3 hours, 15% minutes and 42.5 seconds in the afternoon in the time zone% which is six hours to the west of Coordinated Universal Time% (i.e. Mountain Daylight Time). The time zone specification can% also be written without a colon using one or two-digits% (indicating hours) or three or four digits (indicating hours% and minutes).  Instead of 'seconds' the string may contain 'minutes',% 'hours', 'days' and 'weeks' and all of these may be singular or plural% and have capital first letters.  The letters 'UTC' or 'UT' are allowed at% the end of the string, but these are ignored.%% The parsing of the units attribute is done by the function PARSETNC. This% function can be called by itself to test that a string is sytactically% correct.%% It is possible to have many different types of calendars but timenc only% implements five at present. These are necessary because there is some% confusion with dates before October 15 1582 when the Gregorian calendar was% introduced. A problem also arises when the reference date in the units% attribute is before this. timenc deals with this by recognising some of the% CF conventions and returns different answers depending on the value of the% calendar attribute of the time-like variable. Also, some numerical models% like to pretend that every year has the same number of days - 365, 366 and% 360 are all used.%    1) calendar == 'standard', 'gregorian' (this is the default).% In this case the relevant calculations are done using the functions% get_calendar_date and get_julian_day which know about both the Julian and% Gregorian calendars and so work back to julian day 0 in the year% -4712. This means that the day after 4 October 1582 is 15 October 1582 as% decreed by Pope Gregory XIII. This is the calendar almost universally used% today and what udunits works with today. timenc has worked this way since% revision 1.10 in 2000.%    2) calendar == 'proleptic_gregorian'.% In this case the relevant calculations are done using the matlab functions% datenum and datevec which simply extend the way our present calendar works% backwards into the past. This follows international standard ISO 8601. This% means that dates are continuous, i.e., the day after 4 October 1582 is 5% October 1582, which does NOT correspond to historical time anywhere. As well% there is a year zero. timenc used to work this way before revision 1.10 in the% year 2000 and I believe that udunits also did at some stage in the past.%    3) calendar == 'noleap', '365_day'% Here it is assumed that every year has 365 days.%    4) calendar == 'all_leap', '366_day'% Here it is assumed that every year has 366 days.%    5) calendar == '360_day'% Here it is assumed that every year has 360 days and every month has 30 days.%% Note that other values of the calendar attribute produce an error% message. This can usually be avoided by the user specifying the calendar% explicitly in the call to timenc.%% INPUT:%  file may be the name of a netCDF file with or without the .cdf or .nc%    extent. file may also be the URL of a DODS/OPEnDAP dataset.% time_var: the name of the 'time' variable in the netCDF file or DODS/OPEnDAP%    dataset.  If this argument is missing then it is assumed that variable%    name is 'time'. If time_var is multi-dimensional then it will be handled%    as if it had been reshaped as one 'giant' vector.% bl_corner: a vector of length n specifying the hyperslab corner%    with the lowest index values (the bottom left-hand corner in a%    2-space).  The corners refer to the dimensions in the same%    order that these dimensions are listed in the relevant questions%    in getnc.m and in the inqnc.m description of the variable.  A%    negative element means that all values in that direction will be%    returned.  If this argument is missing or a negative scalar or empty%    array is used this means that all of the elements in the array will be%    returned.%  tr_corner is a vector of length n specifying the hyperslab corner%    with the highest index values (the top right-hand corner in a%    2-space).  The corners refer to the dimensions in the same order%    that these dimensions are listed in the relevant questions in%    getnc.m and in the inqnc.m description of the variable. A negative%    element means that the returned hyperslab should run to the highest%    possible index. Note, however, that the value of an element in the%    tr_corner vector will be ignored if the corresponding element in the%    corner vector is negative.%  calendar is a string determining the type of calendar to be used and is%    discussed above.%% OUTPUT:% gregorian_time: an Mx6 matrix where the rows refer to the M times%    specified in the 'time' variable in the netCDF file.  The columns%    are the year, month, day, hour, minute, second in that order, UT.% serial_time: an M vector giving the serial times (in UT) specified in the%    'time' variable in the netCDF file. Serial times are used by datestr,%    datevec & datenum. Thus gregorian_time = datevec(serial_time). Note that%    the 'time' variable actually contains the serial time relative to a%    base time.% gregorian_base: a 6-vector giving the year, month, day, hour, minute,%    second of the base time as specified in the 'units' attribute of%    the 'time' variable. This is in UT.% serial_base: the serial time of the base time, in UT, as determined by%    matlab's datenum function. Thus gregorian_base = datevec(serial_base).%    serial_base will be a NaN for times before October 15 1582, when the%    Gregorian calendar was adopted, since datenum is not meaningful in this%    case.% sizem: the size of the 'time' variable in the netCDF file.% serial_time_jd: an M vector giving the julian day number (in UT) specified%    in the 'time' variable in the netCDF file. (julian day numbers are used%    by get_julian_day and get_calendar_date. Thus gregorian_time =%    get_calendar_date(serial_time_jd).% serial_base_jd: the Julian day number of the base time, in UT, as%    determined by get_julian_day. Thus gregorian_base =%    get_calendar_date(serial_base_jd).%% EXAMPLES:%       1)%   [gregorian_time, serial_time] = timenc('file.nc', 't');%% This gives you the time information in gregorian and 'matlab' serial% time for a netcdf file named file.nc containing a time variable t.%%       2)%   [gregorian_time, serial_time] = timenc('file.nc', 't', 3, 8);%% This gives you the time information in gregorian and 'matlab' serial% time for elements 3 to 8 of the time variable t.%%       3)%   [gregorian_time, serial_time] = timenc('file.nc', 't', -1, -1, 'gregorian');%% The use of the 5th (calendar) input allows you to override the calendar% attribute in the original netcdf file. This attribute may be invalid or% simply not supported by timenc.%%     Copyright J. V. Mansbridge, CSIRO, Tue May  9 11:36:06 EST 1995%     (with the ability to handle a multi-dimensional time variable%     added by Rose O'Connor).% This function calls: check_nc.m, get_calendar_date.m, loaddap or%                      loaddods, mexnc, parsetnc.m, pos_cds.m%$Id: timenc.m Mon, Wed Nov 29 16:35:09 EDT 2006 $% Process the input argumentsif ~ismember(nargin , [1 2 4 5])  error('timenc takes either 1, 2, 4 or 5 input arguments (try "help timenc")')endif nargin == 1  time_var = 'time';endif nargin <= 2  get_all = 1;elseif nargin >= 4  if isempty(bl_corner) | isempty(tr_corner)    get_all = 1;  else    if (bl_corner < 1)      get_all = 1;    else      get_all = 0;    end  endendif nargin < 5  calendar = [];end% Do some initialisation.CSIRO_add_jar_file_maybe;[mex_name, full_name, desc_das, file_status, exe_name] = ...    choose_mexnc_opendap(file);% We also use method_of_call == 2 which, according to the documentation, does% the same thing as method_of_call == 1 but would be expected to be% slower. However, this is necessary for loaddods to work properly on Windows% boxes when an array is filled with characters. I have no idea what the bug% might be in loaddods. There is also a hint that the same problem can occur% with loaddap on some linux boxes.  method_of_call = 2;switch mex_name case 'mexnc'    %Open the netcdf file.  [cdfid, rcode] = mexnc('OPEN', full_name, 'NC_NOWRITE');  if rcode == -1    error(['** ERROR ** ncopen: rcode = ' num2str(rcode)])  end  %Suppress all error messages from netCDF   [rcode] = mexnc('setopts', 0);  %Get the id number of the variable 'time' and find out the info. about  %the variable 'time'.  [varid, rcode] = mexnc('varid', cdfid, time_var);  if rcode == -1    error(['** ERROR ** ncvarid: time variable = ''' time_var ''' not found'])  end  if get_all == 1    [varnam, vartypv, nvdims, vdims, nvatts, rcode] = ...	mexnc('ncvarinq', cdfid, varid);    if rcode == -1      error(['** ERROR ** ncvarinq: rcode = ' num2str(rcode) ...	     ', time variable = ''' time_var ''''])    end        %Find out the size of the dimension 'time' which is the same as the    %variable 'time'.    %[name, sizem, rcode] = mexnc('ncdiminq', cdfid, vdims(1));    %ROC    sizem = zeros(1,nvdims);    st    = zeros(1,nvdims);    for i = 1:nvdims      [name, sizem(i), rcode] = mexnc('ncdiminq', cdfid, vdims(i));      if rcode == -1	error(['** ERROR ** ncdiminq: rcode = ' num2str(rcode)])      end    end        %Retrieve the elements of the variable 'time'.  These are the serial day    %relative to the base date.    if any(sizem == 0)      serial_rel = []; % Presumably time is an unlimited dimension of length 0.      disp(['Warning: There are apparently no ''time'' records'])    else      %[serial_rel, rcode] = mexnc('ncvarget', cdfid, varid, 0, sizem);      [serial_rel, rcode] = mexnc('ncvarget', cdfid, varid, st, sizem);      if rcode == -1	error(['** ERROR ** ncvarget: rcode = ' num2str(rcode)])      end    end  else    [varnam, vartypv, nvdims, vdims, nvatts, rcode] = ...	mexnc('ncvarinq', cdfid, varid);    if rcode == -1      error(['** ERROR ** ncvarinq: rcode = ' num2str(rcode) ...	     ', time variable = ''' time_var ''''])    end        sizem = zeros(1,nvdims);    for i = 1:nvdims      [name, sizem(i), rcode] = mexnc('ncdiminq', cdfid, vdims(i));      if rcode == -1	error(['** ERROR ** ncdiminq: rcode = ' num2str(rcode)])      end      if tr_corner(i) < 0	tr_corner(i) = sizem(i);      end    end        if any(sizem == 0)      serial_rel = []; % Presumably time is an unlimited dimension of length 0.      disp(['Warning: There are apparently no ''time'' records'])    elseif any(bl_corner < 1)      error('bl_corner values are too small')    elseif any(tr_corner > sizem)      error('tr_corner values are too large')    else      [serial_rel, rcode] = mexnc('ncvarget', cdfid, varid, ...				  bl_corner-1, tr_corner-bl_corner+1);      if rcode == -1	error(['** ERROR ** ncvarget: rcode = ' num2str(rcode)])      end    end  end  %Get the string describing the base date.  [base_str, rcode] = mexnc('attget', cdfid, varid, 'units');  if rcode == -1    error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])  end  %Close the netcdf file.  [rcode] = mexnc('ncclose', cdfid);  if rcode == -1    error(['** ERROR ** ncclose: rcode = ' num2str(rcode)])  end   case {'loaddap', 'loaddods'}  % Dealing with a DODS file.    % Check the value of time_var.   if ischar(time_var)    if ~isfield(desc_das, time_var)      error([time_var ' is not a variable in ' full_name])    end  else    error('For dods data varid must be a string')  end    % Get the structure of time_var and then use this to extract required  % information about attributes. We also strip off extraneous quote marks.    %desc_varid = desc_das.(time_var);  desc_varid = getfield(desc_das, time_var);  if isfield(desc_varid, 'units')    %base_str = desc_varid.('units');    base_str = getfield(desc_varid, 'units');    if ischar(base_str)      base_str = base_str(2:(end - 1));    else      error(['The units attribute of ' time_var ' is not a string'])    end  else    error([time_var ' does not have a units attribute']);  end  if isfield(desc_varid, 'DODS_ML_Size')    %sizem = desc_varid.('DODS_ML_Size');    sizem = getfield(desc_varid, 'DODS_ML_Size');  else    %xx = desc_varid.(time_var);    xx = getfield(desc_varid, time_var);    %sizem = xx.('DODS_ML_Size');    sizem = getfield(xx, 'DODS_ML_Size');  end  num_dim = length(sizem);  % Get the time vector. Note that we allow for the vector of variable values  % to be at different "depths" in the structure. It seems that a variable  % that is also a dimension will be at values_struct.time but otherwise it  % will be at values_struct.time.time.    if get_all    switch method_of_call     case 1      switch mex_name       case 'loaddap'	values_struct = loaddap('+v', [full_name '?' time_var]);       case 'loaddods'	values_struct = loaddods('+v', [full_name '?' time_var]);      end     case 2      switch mex_name       case 'loaddap'	loaddap('+v', [full_name '?' time_var]);       case 'loaddods'	loaddods('+v', [full_name '?' time_var]);      end      eval(['values_struct = ' time_var ';'])    end  else    str_h = [];    for ii = 1:num_dim