function nc_write( outfname, record, fdat, nam, i4, r4, history, long_name )

% function nc_write( outfname, record, fdat, nam, i4, r4, history, long_name );
%
% Creates (if necessary) and writes one record to a netcdf file.
% The inputs are written to mimick the "write_pkd" function.
%
% Return -1 on error, 0 on success
%
% outfname	 output file name, in single quotes
% record	 record number
% fdat		 the data, written as a vector
% nam            character string to describe the variable
% i4, r4         vectors containing metadata information 
%                (same as i4 and r4 vectors used by Alexis Lau)
% history        (optional) Label for additional information (data source,
%                your name, ...).
%                The format is 'string' with single quotes, and it can
%                be as long as you like.  You can either not worry about
%                carriage returns or you can include "char(10)" (no quotes)
%                to force there to be carriage returns.
% long_name      (optional) Single line title for data set.  This is used
%                both as a long_name for the variable and as the title,
%                which is a global attribute.
%
% i4 and r4 comprise 30 and 20 elements, respectively.
% The following elements of i4 and r4 need to be defined:
% i4(1)   year		( 0 for climatology )  
%   (2)   month		
%   (3)   day		
%   (4)   hour		
%   (5)   minutes caution: Alexis" format does not provide for minutes or 
%                 seconds
%   (6)   number of longitude points along a latitude  circle
%   (7)   number of latitude  points along a longitude circle
% r4(1:2) add_offset and scale_factor **** see below ****
%   ( 9)  first longitude point on a latitude  circle
%   (10)  first latitude  point on a longitude circle
%   (11)  longitude resolution
%   (12)  latitude  resolution
%
% 22 October 2013: Modified to run on Matlab R2013b (8.2.0.701)
%
% ***** There is some tricky stuff with add_offset and scale_factor that
% the user needs to be aware of.  add_offset and scale_factor are used to
% pack the data in the following manner.
%	packed = ( unpacked - add_offset ) / scale_factor .
% The values for add_offset and scale_factor are specified are input in r4
% or can be calculated by the algorithm, and are the same for all maps.
% This practice is different than for Alexis Lau's "pkd" format, for which
% add_offset and scale_factor values are determined for each map.
% Special case: for add_offset = scale_factor = NaN, the algorithm calculates
% add_offset and scale_factor from the values of the first map.
%
% Time is specified as "so many hours since some specified time" where
% hours may be fractions of hours.  This choice should placate the 
% GrADS "sdfopen" command and is consistent with the COARDS standard.
%
% Orginal code by Christian Bantzer.
% Modifications by Todd Mitchell, December 1995.
% Time variable made COARDS compliant, October 1996, Todd Mitchell.
% History variable added in October 1999.  Todd Mitchell
% Time units defined to be "hours" (GrADS "sdfopen").  Todd Mitchell 2/2000
% Changed time converstions from mexeps to my juliantime.m 3/2011

variable_name = nam;
if nargin==8; variable_name = long_name; end
variable_units = nam;

xgrid = r4(9)  + [0:i4(6)-1] * r4(11);
ygrid = r4(10) + [0:i4(7)-1] * r4(12);

% If not a row vector, convert to a row vector.
if size(fdat,1)~=1 
[ nrow ncol ] = size( fdat );
fdat = reshape( fdat', 1, nrow*ncol );
end

% The first and subsequent records are handled differently.
% The first record: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if record == 1

% Write time "units" attribute. ---------------------------------

% An older practice is to set the year and perhaps other time variables
% to -999 for a climatology.  Set all time variables < 0 to zero.
x = find(i4(1:5)<0);
if ~isempty(x); i4(x) = zeros(size(x)); end

% The month cannot be zero.
if i4(2)==0 ; 
  disp( sprintf( 'The month cannot be zero.  The month is now 1.' ) );
  i4(2) = 1 ;
end

% The units of time will be specified to be hours so that the GrADS
% "sdfopen" utility will work.
% April 2010: set the reference year to be 1 January 1800.
time_units_string = 'hours since 1800-1-1 0:0:0';
if i4(1)==1;   % A climatology or file where the year doesn't matter.
time_units_string = sprintf( 'hours since %i-%i-%i %i:%i:0', i4(1:5) );
end
% By this convention, the first record has a time value of zero.
% date = 0;   % *** consistent with setting the reference time to 1800.

minutes = 0;
if i4(5)>0 ; 
minutes = fix(i4(5)/60);
i4(5)   = i4(5) - minutes * 60;
end
realtime = i4(1) * 1.0e10 + i4(2) * 1.0e8 + i4(3) * 1.0e6 ...
         + i4(4) * 1.0e4  + minutes*1.0e2 + fix(i4(5));
basetime = 1800  * 1.0e10 + 1.0e8 + 1.0e6;

% Replace the mexeps statement with Matlab functions.
bb = i4(1:5)' ;
time_units_string; 
date = juliantime( bb, time_units_string );
% disp( sprintf( 'i4(1:5) and date are;\n %4d %4d %4d %4d %4d %10d', ...
%   bb, date ) );
% date = mexeps( 'totarray19', realtime, 'hours', basetime );

% Open the netCDF file
cmode = netcdf.getConstant('NETCDF4')
cmode = bitor( cmode, netcdf.getConstant( 'CLASSIC_MODEL' ) )
cdfid = netcdf.create( outfname, cmode )
[ 3 cdfid ]
% cdfid = netcdf.create( outfname, netcdf.getConstant('NETCDF4') ); % this works !!! 

% cdfid = netcdf.create( outfname, 'CLOBBER'); % Really old version
if (cdfid < 0); nc_write = -1; return; end

% Define variable dimensions and variables to hold dimension data
latdimid  = netcdf.defDim( cdfid, 'lat', i4(7) );
londimid  = netcdf.defDim( cdfid, 'lon', i4(6) );
timdimid  = netcdf.defDim( cdfid, 'time', netcdf.getConstant('UNLIMITED') );
latid     = netcdf.defVar( cdfid, 'lat', 'FLOAT', latdimid );
lonid     = netcdf.defVar( cdfid, 'lon', 'FLOAT', londimid );
timid     = netcdf.defVar( cdfid, 'time', 'DOUBLE', timdimid );
% Define a second time variable: yyyymmddhhmm.  *** I was concerned
% about how this would interpreted by other software, so I commented
% it out. T Mitchell, September 2002 ***
% tim2id   = netcdf.defVar( cdfid, 'human_time', 'DOUBLE', 1, timdimid );

% Define the variable.
varid = netcdf.defVar( cdfid, 'data', 'SHORT', ...
  [londimid latdimid timdimid ] );
% varid = netcdf.defVar( cdfid, 'data', 'SHORT', 3, ...
%  [timdimid latdimid londimid]);  % The old ordering. Pay attention ***

if ( (latdimid < 0) | (londimid < 0) ...
  | (timdimid < 0)  | (latid < 0)    | (lonid < 0) ...
  | (timid < 0)     | (varid < 0) )
	nc_write = -1;
	return;
end

% Input the values of add_offset and scale_factor.  If necessary, calculate
% values from the first map.
add_offset   = r4(1);
scale_factor = r4(2);
if isnan(add_offset) & isnan(scale_factor)
x = find( ~isnan( fdat(1,:) ) );
if isempty(x)
disp( sprintf( '\n\nThe first map contains only "NaN"s.' ) );
error( 'You need to input the add_offset and scale_factor (in r4).' );
end
% The scale_factor is chosen to make the scaled number retain 4 significant
% digits (subject to rounding error).
maxval = max( abs( fdat(1,x)-ones(1,length(x))*add_offset ) );
degree = ceil( real( log10( maxval ) ) );
scale_factor = exp( ( degree - 4 ) * log(10) );
disp( sprintf( ...
  '\n\n The add_offset and scale_factor are calculated to be %d %d', ...
  add_offset, scale_factor ) );
end

%Define attributes
% The "Variable not found" error happens in the next line.
netcdf.putAtt( cdfid, varid, 'add_offset',    add_offset );
netcdf.putAtt( cdfid, varid, 'scale_factor',  scale_factor );
% netcdf.putAtt( cdfid, varid, 'missing_value', '32767' ); 
netcdf.defVarFill( cdfid, varid, false,   32767 );
netcdf.putAtt( cdfid, varid, 'units',     variable_units );
netcdf.putAtt( cdfid, varid, 'long_name', variable_name );

netcdf.putAtt( cdfid, timid, 'title', 'Time');
netcdf.putAtt( cdfid, timid, 'units', time_units_string );
netcdf.putAtt( cdfid, timid, 'scale_factor', 1.0 );
netcdf.putAtt( cdfid, timid, 'add_offset',   0.0 );

% netcdf.putAtt( cdfid, tim2id, 'title', 'Human_Time');
% netcdf.putAtt( cdfid, tim2id, 'units', 'yyyymmddhhmm' );

netcdf.putAtt( cdfid, latid, 'title', 'Latitude');
netcdf.putAtt( cdfid, latid, 'units', 'degrees_north');
netcdf.putAtt( cdfid, latid, 'scale_factor', 1.0 );
netcdf.putAtt( cdfid, latid, 'add_offset',   0.0 );

netcdf.putAtt( cdfid, lonid, 'title', 'Longitude');
netcdf.putAtt( cdfid, lonid, 'units', 'degrees_east');
netcdf.putAtt( cdfid, lonid, 'scale_factor', 1.0 );
netcdf.putAtt( cdfid, lonid, 'add_offset',   0.0 );

% Define the "history" label: source of field, netCDF standard, 
% file author, data first record written.
history2 = [ ...
  'The file was written at the JISAO at the University of Washington.', ...
  char(10), ...
'http://jisao.washington.edu' ];
if nargin>=7; history = [ char(10), history, char(10), history2 ]; 
else
history = history2;
end
netcdf.putAtt( cdfid, netcdf.getConstant('NC_GLOBAL'), 'history', history );
netcdf.putAtt( cdfid, netcdf.getConstant('NC_GLOBAL'), 'title',   variable_name );
convention = [ 'The file is written in COARDS-compliant netCDF:', char(10), ...
'ftp://ftp.unidata.ucar.edu/pub/netcdf/Conventions/COARDS/coards_conventions' ];
netcdf.putAtt( cdfid, netcdf.getConstant('NC_GLOBAL'), 'convention', ...
  convention );

% End definition mode
netcdf.endDef( cdfid );

% Write out the latitude and longitude data.
netcdf.putVar( cdfid, latid, ygrid);
netcdf.putVar( cdfid, lonid, xgrid);

netcdf.sync( cdfid);

% Write out the actual data
if (isempty(fdat)) return; end

% Pack the data.
fdat = ( fdat - add_offset ) / scale_factor;

% The next line replaces the NaNs with the value of 32767.  
fdat(isnan(fdat)) = 32767; 
% fdat = fill_nan(fdat, 32767 );

netcdf.putVar( cdfid, timid, [record-1], date );
% netcdf.putVar( cdfid, tim2id, [record-1], ...
%   i4(1)*1.0e8+i4(2)*1.0e6+i4(3)*1.0e4+i4(4)*1.0e2+i4(5) );
netcdf.putVar( cdfid, varid, [0 0 record-1 ], [i4(6) i4(7) 1], fdat);
% netcdf.putVar( cdfid, varid, [record-1 0 0], [1 i4(7) i4(6)], fdat); % Old convention

netcdf.sync( cdfid);

% Subsequent records: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
% The following is for appending a map to an existing netcdf file.
cdfid = netcdf.open( outfname, 'WRITE' );

timid  = netcdf.inqVarID( cdfid, 'time' );
% tim2id = netcdf.inqVarID( cdfid, 'human_time' );
varid  = netcdf.inqVarID( cdfid, 'data' );

time_units_string = netcdf.getAtt( cdfid, timid, 'units' );
if     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),'day') ; time_units = 'days'    ;
else   error( 'The time units is not seconds, hours, or days.' ); end

% test:   time_units_string = 'hours   since    1- 1- 0 00:00:00'
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;
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' );

% Write the input time value for this map in yyyymmddhhmmss form ------

%disp( sprintf( 'i4(1:5): %d %d %d %d %d', i4(1:5) ) );
% An older practice is to set the year and perhaps other time variables
% to -999 for a climatology.  Set all time variables < 0 to zero.
x = find(i4(1:5)<0);
if ~isempty(x); i4(x) = zeros(size(x)); end
%if x~=[] ; i4(x) = zeros(size(x)); end

minutes = 0;
if i4(5)>0 ; 
minutes = fix(i4(5)/60);
i4(5)   = i4(5) - minutes * 60;
end

% The month and day cannot be zero.
if i4(3)==0 ; i4(3) = 1 ; end
if i4(2)==0 ; i4(2) = 1 ; end

% realtime = i4(1) * 1.0e10 + i4(2) * 1.0e8 + i4(3) * 1.0e6 ...
%          + i4(4) * 1.0e4  + minutes*1.0e2 + fix(i4(5));
% date = mexeps( 'totarray19', realtime, time_units, basetime );
bb = i4(1:5)';
time_units_string;
date = juliantime( bb, time_units_string );

add_offset   = netcdf.getAtt( cdfid, varid, 'add_offset'   );
scale_factor = netcdf.getAtt( cdfid, varid, 'scale_factor' );
%disp( sprintf('add_offset, scale_factor: %d %d\n', ...
%  add_offset, scale_factor ) );
% Pack the data.
fdat = ( fdat - add_offset ) / scale_factor;
fdat(isnan(fdat)) = 32767; 
% fdat = fill_nan(fdat, 32767 );
%disp( sprintf( 'fdat(1:5) %d %d %d %d %d\n', fdat(1:5) ) )

netcdf.putVar( cdfid, timid, [record-1], 1, date );
% netcdf.putVar( cdfid, tim2id, [record-1], 1, ...
%   i4(1)*1.0e8+i4(2)*1.0e6+i4(3)*1.0e4+i4(4)*1.0e2+i4(5) );
netcdf.putVar( cdfid, varid, [ 0 0 record-1 ], [i4(6) i4(7) 1 ], fdat);
% netcdf.putVar( cdfid, varid, [record-1 0 0], [1 i4(7) i4(6) ], fdat); % Old ordering
netcdf.sync( cdfid);

end

netcdf.close( cdfid);

% disp( sprintf( 'Have successfully written out map %5d of data.', record ) )
nc_write = 0;

return;

	




