%RCMRAW2CDF convert 5059 reader program ascci output to netcdf
% function rcmraw2cdf(rcmfile, ncfile, recnums, aa, metadata);
% rcmfile = raw aanderaa data exported from 5059 reader program
% How to export the data this file reads:
%   do not apply a station, use the pristine raw data right from the DSU.
%   make sure time is displayed in the raw data list
%   under tooling -> wizards, hit the export to ascii wizard
%   using custom settings, tabs between colums, cr-lf between samples, 
%   include time tags and record numbers, include channel order in header
% ncfile = netCDF output file name
% recnums = [start end] samples to trim data to (to eliminate crappy data)
%   set end to inf to get all the data
% aa.serial = RCM serial number
% aa.instheight = mab of the instrument
% aa.nchannels = number of channels recorded
% aa.interval = sample interval in seconds
% aa.chan2.serial = DCS serial number, use '' if there is no DCS
% aa.chan2.cals = [A B C D] = [0 2.933E-1 0 0] for speed
% aa.chan3.cals  = [A B C D] = [0 3.516E-1 0 0] for direction
% aa.chan4.serial = temp serial number, use '' if there is no thermistor
% aa.chan4.range = Wide | Low | High | Arctic
% aa.chan4.cals  = [A B C D] for temp
% aa.chan5.serial = cond serial number, use '' if there is no cond
% aa.chan5.range = 0-74 | 24-38 | 0-2
% aa.chan5.cals  = [A B C D] for cond
% aa.chan6.serial = press serial number, use '' if there is no press
% aa.chan6.cals  = [A B C D] for press
% aa.chan7.serial = turb serial number, use '' if there is no turbidity
% aa.chan7.cals  = [A B C D] for turbidity
% aa.chan8.serial = optode serial number, use '' if there is no optode
% aa.chan8.units = uM/l | %
% aa.chan8.cals  = [A B C D] for optode (or other sensor)
% this program applies the calibrations, user be prompted for the calibration constants if missing
% the metadata
% metadata.mooring_number = '7281'; %four digit NNNL, NNN mooring + L logger
% metadata.deployment_start = '08-nov-2002 12:33:00'; % the in water time
% metadata.deployment_end = '16-feb-2003 13:05:00'; % the out of water time
% metadata.lon = 13.29722; % decimal degrees, West = negative
% metadata.lat = 43.75652; % decimal degrees, South = negative
% metadata.declination = 2.0; % degrees west is negative
% metadata.water_depth = 10.0; % m 
% metadata.metafile = [pwd '\meta7023'];
% metadata.metafile_data = datestr(now,1);
% metadata.metafile_version = '1.0';
% metadata.metafile_author = 'RMH';
% % if the following fields are missing, UNKOWN will be used
% metadata.origin = 'USGS/WHFC'; % with collaborator's data, could be USC, etc.
% metadata.experiment = 'Boston Longterm'; 
% metadata.project = 'WHFC'; % might also use OFA funding agency, such as
% MWRA, EPA, WCMG
% metadata.adjust_time = 1; % use time stamps and interval to recompute
% time
%
% time_stamp = the matrix of daily time_stamps found in the file
%
% Version 2.2 11/11/04
% trap times that are NaNs and use surrounding times to fix them
% Version 2.1 10/29/04
% trying to deal with missing records
% Version 2.0 10/27/04
% remove batch functionality, use metadata structure
% Version 1.0
% add batch functionality for convenient processing with all the coefficients

% Written by Marinna Martini for the Woods Hole Field Center
% 5/24/04   The first real processing tool version

% TODO might need to trap time stamps better
% TODO test with pressure, cond and turbidity sensors

function time_stamp = rcmraw2cdf(rcmfile, ncfile, recnums, aa, metadata)

mver = '2.1';
maxchannels = 8;

% -------- verify inputs
if ~exist('rcmfile','var'),
    [rcmfile, path, filtidx] = uigetfile('*.Asc','Open aanderaa 5059 ASCII export of RCM raw data');
    if isempty(rcmfile), return; end
    rcmfile= [path rcmfile];
end
if ~exist('ncfile','var'), % if not specified, make it up
    [path,name,ext,ver] = fileparts(rcmfile);
    ncfile = fullfile(path, [name '.cdf']);
end
if ~exist('aa','var'), % if not specified, abort
    disp('Aanderaa information and calibrations missing')
    return
end
if ~exist('metadata','var'), % if not specified, abort
    disp('metadata missing')
    return
end

% get all the data at once, these are not huge files
fid = fopen(rcmfile);
disp('Reading the rcm file...')
tic
% skip header
fgetl(fid);
fgetl(fid);

buf = fgetl(fid);
irec = 1;
istamp = 1;
newrec = 0;
stamptime = [];
newstamp = 0;
%for i=1:40, 
while ischar(buf),
    [aarec, count, errmsg] = sscanf(buf,'%d %d/%d/%d %d:%d:%d %1c%*1c %4d %4d %4d %4d %4d %4d %4d %4d');
    if (count == 14) && (aarec(9) == 7), % it's a time stamp...
        newstamp = 1;
        time_stamp(istamp,:) = aarec';
        %disp(sprintf('Time stamp at RCM rec #%d: %s',aarec(1), buf))
        % do the time from the time stamps
        [tm, newstamptime, dt] = aa2jul(time_stamp(istamp,2:8));
        istamp = istamp+1;
    else
        if count == 8+aa.nchannels, % probably a good data record..
            raw(irec,:) = aarec';
        else % it's crap
            disp(sprintf('Bad data at RCM rec #%d: %s',aarec(1), buf))
            % we may want to insert time here
            raw(irec,:) = ones(1,8+aa.nchannels).*NaN;
        end
        if newstamp, stamptime(irec) = newstamptime;
        else stamptime(irec) = stamptime(irec-1)+aa.interval./(24*3600);
        end
        %if irec>1000 && irec<1200, disp(sprintf('%d %s %s',irec, datestr(gregorian(stamptime(irec))), buf)); end
        if irec<1000 && ~rem(irec,100), disp(sprintf('%d samples read in %f min',irec, toc/60)), end
        if irec>1000 && ~rem(irec,1000), disp(sprintf('%d samples read in %f min',irec, toc/60)), end
        irec = irec + 1;
        newstamp = 0;
    end
    buf = fgetl(fid);
end
fclose(fid);

% trim the bad data
if ~exist('recnums','var'), recnums = [1 length(raw)]; end
if isinf(recnums(1)) || (recnums(1) < 1), recnums(1) = 1; end
if isinf(recnums(2)) || (recnums(2) > length(raw)), recnums(2) = length(raw); end
% find the record number in the file records
irec = find(raw(:,1) == recnums(1));
if isempty(irec), irec = 1; 
elseif length(irec) > 1, irec = min(irec);
end
recnums(1) = irec;
irec = find(raw(:,1) == recnums(2));
if isempty(irec), irec = length(raw(:,1)); 
elseif length(irec) > 1, irec = min(irec);
end
recnums(2) = irec;

raw = raw(recnums(1):recnums(2),:);
stamptime = stamptime(recnums(1):recnums(2));

% apply the calibrations
if aa.nchannels > 1,
    aa.chan2.data = aa.chan2.cals(1)+aa.chan2.cals(2).*raw(:,10)+aa.chan2.cals(3).*(raw(:,10).^2)+aa.chan2.cals(4).*(raw(:,10).^3);
end
if aa.nchannels > 2,
    aa.chan3.data = aa.chan3.cals(1)+aa.chan3.cals(2).*raw(:,11)+aa.chan3.cals(3).*(raw(:,11).^2)+aa.chan3.cals(4).*(raw(:,11).^3);
end
if aa.nchannels > 3,
    aa.chan4.data = aa.chan4.cals(1)+aa.chan4.cals(2).*raw(:,12)+aa.chan4.cals(3).*(raw(:,12).^2)+aa.chan4.cals(4).*(raw(:,12).^3);
end
if aa.nchannels > 4,
    aa.chan5.data = aa.chan5.cals(1)+aa.chan5.cals(2).*raw(:,13)+aa.chan5.cals(3).*(raw(:,13).^2)+aa.chan5.cals(4).*(raw(:,13).^3);
end
if aa.nchannels > 5,
    aa.chan6.data = aa.chan6.cals(1)+aa.chan6.cals(2).*raw(:,14)+aa.chan6.cals(3).*(raw(:,14).^2)+aa.chan6.cals(4).*(raw(:,14).^3);
end
if aa.nchannels > 6,
    aa.chan7.data = aa.chan7.cals(1)+aa.chan7.cals(2).*raw(:,15)+aa.chan7.cals(3).*(raw(:,15).^2)+aa.chan7.cals(4).*(raw(:,15).^3);
end
if aa.nchannels > 7,
    aa.chan8.data = aa.chan8.cals(1)+aa.chan8.cals(2).*raw(:,16)+aa.chan8.cals(3).*(raw(:,16).^2)+aa.chan8.cals(4).*(raw(:,16).^3);
end

cdf = netcdf(ncfile, 'clobber');
if isempty(cdf), return, end

definencfile(cdf,aa,metadata);
cdf.INST_TYPE = ncchar('RCM-9MkII');

add_vardesc(cdf);
add_fillvalues(cdf,1E35);

endef(cdf);

% write the variables.  Note... must do this after endef!
if metadata.adjust_time, % do the time from the time stamps and interval
    cdf{'time'}(1:length(stamptime)) = floor(stamptime);
    cdf{'time2'}(1:length(stamptime)) = (stamptime-floor(stamptime)).*(1000*24*3600);    
    idx = min(find(~isnan(stamptime)));
    cdf.start_time = ncchar(datestr(datenum(gregoria(stamptime(idx)))));
    idx = max(find(~isnan(stamptime)));
    cdf.stop_time = ncchar(datestr(datenum(gregoria(stamptime(idx)))));
    cdf.DELTA_T = ncchar(aa.interval);
else % do the time from the good data records
    [tm, tj, dt] = aa2jul(raw(:,2:8));
    % find any NaN's and fix them, time only
    tj = fix_time_nans(tj,aa.interval);
    cdf{'time'}(1:length(tj)) = floor(tj);
    cdf{'time2'}(1:length(tj)) = (tj-floor(tj)).*(1000*24*3600);
    cdf.start_time = ncchar(datestr(datenum(gregoria(tj(1)))));
    cdf.stop_time = ncchar(datestr(datenum(gregoria(tj(end)))));
    cdf.DELTA_T = ncchar(floor(int2str(dt)));
end
cdf{'depth'}(1) = metadata.water_depth;
cdf{'lon'}(1) = metadata.lon;
cdf{'lat'}(1) = metadata.lat;
cdf{'recnum'}(1:length(raw(:,1))) = raw(:,1);
if ~isempty(aa.chan2.serial), cdf{'CS_300'}(:) = aa.chan2.data; end
if ~isempty(aa.chan3.serial), 
    dirdat = aa.chan3.data (:) + metadata.declination;
    beyond = find(dirdat<0);
    dirdat(beyond) = dirdat(beyond) + 360;
    beyond = find(dirdat>360);
    dirdat(beyond) = dirdat(beyond) - 360;
    cdf{'CD_310'}(:) = dirdat(:); 
    dir = dirdat(:) .* pi ./180.;
    [v,u] = pol2cart(dir,aa.chan2.data); 
    cdf{'u_1205'}(:) = u;
    cdf{'v_1206'}(:) = v;
end
if ~isempty(aa.chan4.serial), cdf{'T_28'}(1:length(raw(:,1))) = aa.chan4.data; end
if ~isempty(aa.chan5.serial), cdf{'C_51'}(1:length(raw(:,1))) = aa.chan5.data./10; end % convert from mS/cm to S/m
if ~isempty(aa.chan6.serial), cdf{'P_1'}(1:length(raw(:,1))) = aa.chan6.data.*10; end % convert from kPa to db
% 1db = 1x10^4kg/(ms^2) and 1kPa = 1x10^3kg/(ms^2)
if ~isempty(aa.chan7.serial), cdf{'nep'}(:) = aa.chan7.data; end 

if ~isempty(aa.chan8.serial), 
    if strcmp('uM/l',aa.chan8.units), 
        % use ml/l at this stage because we do not yet have salinity to
        % calculate density and therefore uM/kg
        % RCM outputs data in counts, uM/l = counts * 0.488281
        % 0.488281 was multiplied as one of the cals
        % ml/l = uM/l / 44.66
        cdf{'O_60'}(1:length(raw(:,1))) = aa.chan8.data./44.66;
    else cdf{'O_62'}(1:length(raw(:,1))) = aa.chan8.data;
    end
    cdf{'o2raw'}(1:length(raw(:,1))) = raw(:,16);
end

cdf.history = ncchar(sprintf('Created by rcmraw2cdf.m %s, declination of %f applied.',mver,metadata.declination));
add_minmaxvalues(cdf);

disp(sprintf('Wrote %d records from %s to %s', length(cdf{'time'}(:)),...
    cdf.start_time(:), cdf.stop_time(:)))

close(cdf);

disp(sprintf('Rcmraw2cdf done in %f min',toc/60))

return % from main

% ********************************** subfunctions ********************************

function [tm, tj, dt] = aa2jul(aaTime)
% aaTime should be read from the file in the format
% 1 2 3    4 5  6  7
% 2/5/2004 9:11:05 PM
% construct the time from the time stamp string
% first, figure out AMs and PMs
    pmidx = find((aaTime(:,4) ~= 12) & aaTime(:,7) == 80);
    aaTime(pmidx,4) = aaTime(pmidx,4)+12;   % convert afternoon times to 24 hour format
    amidx = find((aaTime(:,4) == 12) & (aaTime(:,7) == 65));  % of the AM times, find the 12 AM instances
    aaTime(amidx,4) = aaTime(amidx,4)-12;   % convert 12:** AM format to 00:** hours
    tm = datenum(aaTime(:,3),aaTime(:,1),aaTime(:,2),aaTime(:,4),aaTime(:,5),aaTime(:,6));
    tj = julian([aaTime(:,3) aaTime(:,1) aaTime(:,2) aaTime(:,4) aaTime(:,5) aaTime(:,6)]);
    if length(tj) > 1, dt = mean(diff(tj)).*(24*3600);
    else dt = []; end
return

function aa = get_meta()
    sensor = 'Aanderaa RCM-9 MkII Information';
    prompt  = {'Enter the RCM serial number:',...
        'Enter height above bottom, m:',...
        'Enter number of channels recorded:',...
        'Enter depth at the site:'};
    def     = {'RCM###','0','1','0'};
    dlgtitle   = ['Input metadata for ' sensor];
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    aa.serial = dlgresult{1};
    aa.instheight = str2num(dlgresult{2});
    aa.nchannels = str2num(dlgresult{3});

    if aa.nchannels > 1,
        sensor = 'Aanderaa DCS 3920S Current Speed';
        units = 'cm/s';
        prompt  = {'Enter the DCS serial number on channel 2 (use '''' for no sensor):',...
            'Enter Coefficient A:',...
            'Enter Coefficient B:',...
            'Enter Coefficient C:',...
            'Enter Coefficient D:'};
        def     = {'','0','2.933e-01','0','0'};
        dlgtitle   = ['Input metadata for ' sensor];
        lineNo  = 1;
        dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
        aa.chan2.serial = dlgresult{1};
        aa.chan2.cals(1) = str2num(dlgresult{2});
        aa.chan2.cals(2) = str2num(dlgresult{3});
        aa.chan2.cals(3) = str2num(dlgresult{4});
        aa.chan2.cals(4) = str2num(dlgresult{5});
    end

    if aa.nchannels > 2,
    sensor = 'Aanderaa DCS 3920S Current Direction';
    units = 'Deg. M';
    prompt  = {'Enter the calibration date (use '''' for no sensor):',...
        'Enter the DCS serial number:',...
        'Enter Coefficient A:',...
        'Enter Coefficient B:',...
        'Enter Coefficient C:',...
        'Enter Coefficient D:'};
    def     = {'','0','3.516E-01','0','0'};
    dlgtitle   = ['Input metadata for ' sensor];
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    aa.chan3.serial = dlgresult{1};
    aa.chan3.cals(1) = str2num(dlgresult{2});
    aa.chan3.cals(2) = str2num(dlgresult{3});
    aa.chan3.cals(3) = str2num(dlgresult{4});
    aa.chan3.cals(4) = str2num(dlgresult{5});
    end

    if aa.nchannels > 3,
    sensor = 'Aanderaa Temperature Model 3621';
    units = 'deg. C';
    prompt  = {'Enter the sensor serial number (use '''' for no sensor):',...
        'Enter the temperature range:',...
        'Enter Coefficient A:',...
        'Enter Coefficient B:',...
        'Enter Coefficient C:',...
        'Enter Coefficient D:'};
    def     = {'1105', 'low','-2.761e00','2.406e-02','-2.238e-06','2.056e-09'};
    dlgtitle   = ['Input metadata for ' sensor];
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    aa.chan4.serial = dlgresult{1};
    aa.chan4.range = dlgresult{2};
    aa.chan4.cals(1) = str2num(dlgresult{3});
    aa.chan4.cals(2) = str2num(dlgresult{4});
    aa.chan4.cals(3) = str2num(dlgresult{5});
    aa.chan4.cals(4) = str2num(dlgresult{6});
    end

    if aa.nchannels > 4,
    sensor = 'Aanderaa Conductivity';
    units = 'mS/cm';
    prompt  = {'Enter the sensor serial number (use '''' for no sensor):',...
        'Enter the conductivity range:',...
        'Enter Coefficient A:',...
        'Enter Coefficient B:',...
        'Enter Coefficient C:',...
        'Enter Coefficient D:'};
    def     = {'', '0-74','0','0','0','0'};
    dlgtitle   = ['Input metadata for ' sensor];
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    aa.chan5.serial = dlgresult{1};
    aa.chan5.range = dlgresult{2};
    aa.chan5.cals(1) = str2num(dlgresult{3});
    aa.chan5.cals(2) = str2num(dlgresult{4});
    aa.chan5.cals(3) = str2num(dlgresult{5});
    aa.chan5.cals(4) = str2num(dlgresult{6});
    end

    if aa.nchannels > 5,
    sensor = 'Aanderaa Pressue';
    units = 'KPa/MPa';
    prompt  = {'Enter the sensor serial number (use '''' for no sensor):',...
        'Enter Coefficient A:',...
        'Enter Coefficient B:',...
        'Enter Coefficient C:',...
        'Enter Coefficient D:'};
    def     = {'','0','0','0','0'};
    dlgtitle   = ['Input metadata for ' sensor];
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    aa.chan6.serial = dlgresult{1};
    aa.chan6.cals(1) = str2num(dlgresult{2});
    aa.chan6.cals(2) = str2num(dlgresult{3});
    aa.chan6.cals(3) = str2num(dlgresult{4});
    aa.chan6.cals(4) = str2num(dlgresult{5});
    end

    if aa.nchannels > 6,
    sensor = 'Aanderaa Turbidity';
    units = 'NTU';
    prompt  = {'Enter the sensor serial number (use '''' for no sensor):',...
        'Enter Coefficient A:',...
        'Enter Coefficient B:',...
        'Enter Coefficient C:',...
        'Enter Coefficient D:'};
    def     = {'','0','0','0','0'};
    dlgtitle   = ['Input metadata for ' sensor];
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    aa.chan7.serial = dlgresult{1};
    aa.chan7.cals(1) = str2num(dlgresult{2});
    aa.chan7.cals(2) = str2num(dlgresult{3});
    aa.chan7.cals(3) = str2num(dlgresult{4});
    aa.chan7.cals(4) = str2num(dlgresult{5});
    end

    if aa.nchannels > 7,
    sensor = 'Aanderaa Optode 3830';
    prompt  = {'Enter the optode serial number (use '''' for no sensor):',...
        'Enter the output units:',...
        'Enter Coefficient A:',...
        'Enter Coefficient B:',...
        'Enter Coefficient C:',...
        'Enter Coefficient D:',...
        'Enter Salinity Setting, ppt:'};
    def     = { '113','uM/l','0','0.488281','0','0','0'};
    dlgtitle   = ['Input metadata for ' sensor];
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    aa.chan8.serial = dlgresult{1};
    aa.chan8.units = dlgresult{2};
    aa.chan8.cals(1) = str2num(dlgresult{3});
    aa.chan8.cals(2) = str2num(dlgresult{4});
    aa.chan8.cals(3) = str2num(dlgresult{5});
    aa.chan8.cals(4) = str2num(dlgresult{6});
    aa.chan8.salinity = str2num(dlgresult{7});
    end
    
    if nchannels < maxchannels,
        for n = nchannels+1:maxchannels, eval(sprintf('aa.chan%d.serial = '''';'),n),        end
    end
return

function definencfile(cdf,aa,metadata)

% EPIC requirements USGS uses
if isfield(metadata,'mooring_number'), cdf.MOORING = ncchar(metadata.mooring_number);
else cdf.MOORING = ncchar('XXXX'); end
if isfield(metadata,'water_depth'), cdf.WATER_DEPTH = ncfloat(metadata.water_depth);
else cdf.WATER_DEPTH = nclong(99); end
if isfield(metadata,'lat'), cdf.latitude = ncfloat(metadata.lat);
else cdf.latitude = ncfloat(999); end
if isfield(metadata,'lon'), cdf.longitude = ncfloat(metadata.lon);
else cdf.longitude = ncfloat(999); end
if isfield(metadata,'declination'), cdf.magnetic_variation = ncfloat(metadata.declination);
else cdf.magnetic_variation = ncfloat(0); end
cdf.INST_TYPE = ncchar('RCM9MkII');
if isfield(metadata,'origin'), cdf.DATA_ORIGIN = ncchar(metadata.origin);
else cdf.DATA_ORIGIN = ncchar('UNKNOWN'); end
if isfield(metadata,'experiment'), cdf.EXPERIMENT = ncchar(metadata.experiment);
else cdf.EXPERIMENT = ncchar('UNKOWN'); end
if isfield(metadata,'project'), cdf.PROJECT = ncchar(metadata.project);
else cdf.PROJECT = ncchar('UNKOWN'); end

% things simply to stay EPIC compliant, not that USGS needs them
cdf.DATA_SUBTYPE = ncchar('');
cdf.WATER_MASS = ncchar('?');
cdf.COMPOSITE = nclong(0); % for flagging patched together data sets
cdf.POS_CONST = nclong(0);
cdf.FILL_FLAG = nclong(1);
cdf.DEPTH_CONST = nclong(1); % depth will change
cdf.DATA_CMNT = ncchar('');
cdf.VAR_FILL = ncfloat(1e35);

cdf.Deployment_date = metadata.deployment_start;
cdf.Recovery_date = metadata.deployment_end;

%% Dimensions:
cdf('time') = 0; %% (record dimension)
cdf('depth') = 1;
cdf('lon') = 1;
cdf('lat') = 1;

%% Variables and attributes:
cdf{'time'} = nclong('time'); 
cdf{'time'}.FORTRAN_format = ncchar('F10.2');
cdf{'time'}.units = ncchar('True Julian Day');
cdf{'time'}.type = ncchar('UNEVEN');
cdf{'time'}.epic_code = nclong(624);
if metadata.adjust_time,
    cdf{'time'}.comment = ncchar('Time calculated from daily time stamps');
else 
    cdf{'time'}.comment = ncchar('Time from record time stamps');
end

cdf{'time2'} = nclong('time'); 
cdf{'time2'}.FORTRAN_format = ncchar('F10.2');
cdf{'time2'}.units = ncchar('msec since 0:00 GMT');
cdf{'time2'}.type = ncchar('EVEN');
cdf{'time2'}.epic_code = nclong(624);

cdf{'depth'} = ncfloat('depth'); %% 1 element.
cdf{'depth'}.FORTRAN_format = ncchar('F10.2');
cdf{'depth'}.units = ncchar('m');
cdf{'depth'}.type = ncchar('EVEN');
cdf{'depth'}.epic_code = nclong(3);
cdf{'depth'}(:) = metadata.water_depth;

cdf{'lon'} = ncfloat('lon'); %% 1 element.
cdf{'lon'}.FORTRAN_format = ncchar('F10.2');
cdf{'lon'}.units = ncchar('degrees_east');
cdf{'lon'}.type = ncchar('EVEN');
cdf{'lon'}.epic_code = nclong(502);
cdf{'lon'}(:) = metadata.lon;
 
cdf{'lat'} = ncfloat('lat'); %% 1 element.
cdf{'lat'}.FORTRAN_format = ncchar('F10.2');
cdf{'lat'}.units = ncchar('degrees_north');
cdf{'lat'}.type = ncchar('EVEN');
cdf{'lat'}.epic_code = nclong(500);
cdf{'lat'}(:) = metadata.lat;

cdf{'recnum'} = ncfloat('time', 'depth', 'lat', 'lon');
cdf{'recnum'}.name = ncchar('recnum');
cdf{'recnum'}.long_name = ncchar('record number          ');
cdf{'recnum'}.generic_name = ncchar('rec');
cdf{'recnum'}.units = ncchar('count');
cdf{'recnum'}.valid_range = ncfloat([0 262100]);
cdf{'recnum'}.serial_number = ncchar(aa.serial);

if ~isempty(aa.chan4.serial),
    cdf{'T_28'} = ncfloat('time', 'depth', 'lat', 'lon');
    cdf{'T_28'}.name = ncchar('T');
    cdf{'T_28'}.long_name = ncchar('TEMPERATURE (C)          ');
    cdf{'T_28'}.generic_name = ncchar('temp');
    cdf{'T_28'}.FORTRAN_format = ncchar('f10.2');
    cdf{'T_28'}.units = ncchar('C');
    cdf{'T_28'}.epic_code = nclong(28);
    cdf{'T_28'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdf{'T_28'}.serial_number = ncchar(aa.chan4.serial);
    cdf{'T_28'}.valid_range = ncfloat([0 100]);
    cdf{'T_28'}.cal_A = ncdouble(aa.chan4.cals(1));
    cdf{'T_28'}.cal_B = ncdouble(aa.chan4.cals(2));
    cdf{'T_28'}.cal_C = ncdouble(aa.chan4.cals(3));
    cdf{'T_28'}.cal_D = ncdouble(aa.chan4.cals(4));
end

if ~isempty(aa.chan2.serial),
    % 300:CS :CURRENT SPEED (CM/S)      :vspd:cm s-1:f8.2:oceanographic (going to)
    cdf{'CS_300'} = ncfloat('time', 'depth', 'lat', 'lon');
    cdf{'CS_300'}.name = ncchar('CS');
    cdf{'CS_300'}.long_name = ncchar('Current Speed');
    cdf{'CS_300'}.generic_name = ncchar('vspd');
    cdf{'CS_300'}.FORTRAN_format = ncchar('f8.2');
    cdf{'CS_300'}.units = ncchar('cm/s');
    cdf{'CS_300'}.epic_code = nclong(300);
    cdf{'CS_300'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdf{'CS_300'}.serial_number = nclong(aa.chan2.serial);
    cdf{'CS_300'}.valid_range = ncfloat([0 1000]);
    cdf{'CS_300'}.cal_A = ncdouble(aa.chan2.cals(1));
    cdf{'CS_300'}.cal_B = ncdouble(aa.chan2.cals(2));
    cdf{'CS_300'}.cal_C = ncdouble(aa.chan2.cals(3));
    cdf{'CS_300'}.cal_D = ncdouble(aa.chan2.cals(4));
end

if ~isempty(aa.chan3.serial),
    % 310:CD :CURRENT DIRECTION (T)     :vdir:degrees:f8.2:oceanographic (going to), using true north
    cdf{'CD_310'} = ncfloat('time', 'depth', 'lat', 'lon');
    cdf{'CD_310'}.name = ncchar('CD');
    cdf{'CD_310'}.long_name = ncchar('Current Direction');
    cdf{'CD_310'}.generic_name = ncchar('vdir');
    cdf{'CD_310'}.FORTRAN_format = ncchar('f8.2');
    cdf{'CD_310'}.units = ncchar('degrees');
    cdf{'CD_310'}.epic_code = nclong(310);
    cdf{'CD_310'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdf{'CD_310'}.serial_number = nclong(aa.chan3.serial);
    cdf{'CD_310'}.valid_range = ncfloat([0 1000]);
    cdf{'CD_310'}.cal_A = ncdouble(aa.chan3.cals(1));
    cdf{'CD_310'}.cal_B = ncdouble(aa.chan3.cals(2));
    cdf{'CD_310'}.cal_C = ncdouble(aa.chan3.cals(3));
    cdf{'CD_310'}.cal_D = ncdouble(aa.chan3.cals(4));
    disp(sprintf('Declination of %f has been added to the current direction to write true heading to the file',...
        metadata.declination));
end

if ~isempty(aa.chan2.serial),
    % 1205:u   :Eastward Velocity        :u:cm/s: :Eastward Velocity
    cdf{'u_1205'} = ncfloat('time', 'depth', 'lat', 'lon'); cdfobj = cdf{'u_1205'};
    cdfobj.name = ncchar('u');
    cdfobj.long_name = ncchar('Eastward Velocity        ');
    cdfobj.generic_name = ncchar('u');
    cdfobj.FORTRAN_format = ncchar('');
    cdfobj.units = ncchar('cm/s');
    cdfobj.epic_code = nclong(1205);
    cdfobj.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdfobj.valid_range = ncfloat([-1000 1000]);
end

if ~isempty(aa.chan3.serial),
    %1206:v   :Northward Velocity       :v:cm/s: :Northward Velocity    
    cdf{'v_1206'} = ncfloat('time', 'depth', 'lat', 'lon'); cdfobj = cdf{'v_1206'};
    cdfobj.name = ncchar('v');
    cdfobj.long_name = ncchar('Northward Velocity       ');
    cdfobj.generic_name = ncchar('v');
    cdfobj.FORTRAN_format = ncchar('');
    cdfobj.units = ncchar('cm/s');
    cdfobj.epic_code = nclong(1206);
    cdfobj.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdfobj.valid_range = ncfloat([-1000 1000]);
end

if ~isempty(aa.chan5.serial),
    % 51:C  :CONDUCTIVITY             :con:S/m:f10.3:
    cdf{'C_51'} = ncfloat('time', 'depth', 'lat', 'lon');
    cdf{'C_51'}.name = ncchar('C');
    cdf{'C_51'}.long_name = ncchar('CONDUCTIVITY');
    cdf{'C_51'}.generic_name = ncchar('con');
    cdf{'C_51'}.FORTRAN_format = ncchar('f10.3');
    cdf{'C_51'}.units = ncchar('S/m');
    cdf{'C_51'}.epic_code = nclong(51);
    cdf{'C_51'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdf{'C_51'}.serial_number = nclong(aa.chan5.serial);
    cdf{'C_51'}.valid_range = ncfloat([0 7.4]);
    cdf{'C_51'}.cal_A = ncdouble(aa.chan5.cals(1));
    cdf{'C_51'}.cal_B = ncdouble(aa.chan5.cals(2));
    cdf{'C_51'}.cal_C = ncdouble(aa.chan5.cals(3));
    cdf{'C_51'}.cal_D = ncdouble(aa.chan5.cals(4));
end

if ~isempty(aa.chan6.serial),
    % 1:P  :PRESSURE (DB)            :depth:dbar:f10.1:
    % 4:P  :PRESSURE (PASCALS)       :depth:Pa: :
    cdf{'P_1'} = ncfloat('time', 'depth', 'lat', 'lon');
    cdf{'P_1'}.name = ncchar('P');
    cdf{'P_1'}.long_name = ncchar('PRESSURE (PASCALS)       ');
    cdf{'P_1'}.generic_name = ncchar('depth');
    cdf{'P_1'}.FORTRAN_format = ncchar('');
    cdf{'P_1'}.units = ncchar('Pa');
    cdf{'P_1'}.epic_code = nclong(4);
    cdf{'P_1'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdf{'P_1'}.serial_number = nclong(aa.chan6.serial);
    cdf{'P_1'}.valid_range = ncfloat([0 20000]);
    cdf{'P_1'}.cal_A = ncdouble(aa.chan6.cals(1));
    cdf{'P_1'}.cal_B = ncdouble(aa.chan6.cals(2));
    cdf{'P_1'}.cal_C = ncdouble(aa.chan6.cals(3));
    cdf{'P_1'}.cal_D = ncdouble(aa.chan6.cals(4));
end

if ~isempty(aa.chan7.serial),
    % probably NTU
    % moddeled after
    % 56:NEP:BACKSCATTER INTENSITY    :nephylometer:v:f10.6: 
    cdf{'nep'} = ncfloat('time', 'depth', 'lat', 'lon');
    cdf{'nep'}.long_name = ncchar('turbidity');
    cdf{'nep'}.units = ncchar('volts');
    cdf{'nep'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdf{'nep'}.serial_number = nclong(aa.chan7.serial);
    cdf{'nep'}.valid_range = ncfloat([0 5]);
    cdf{'nep'}.cal_A = ncdouble(aa.chan7.cals(1));
    cdf{'nep'}.cal_B = ncdouble(aa.chan7.cals(2));
    cdf{'nep'}.cal_C = ncdouble(aa.chan7.cals(3));
    cdf{'nep'}.cal_D = ncdouble(aa.chan7.cals(4));
end

if ~isempty(aa.chan8.serial),
    if strcmp('uM/l',aa.chan8.units),
        % use ml/l at this stage because we do not yet have salinity to
        % calculate density and therefore uM/kg
        % RCM outputs data in counts, uM/l = counts * 0.488281
        % ml/l = uM/l / 44.66
        % 60:O  :OXYGEN (ML/L)            :ox:ml/l:f10.2:Dissolved oxygen calculated from CTD values
        cdf{'O_60'} = ncfloat('time', 'depth', 'lat', 'lon'); cdfobj = cdf{'O_60'};
        cdfobj.name = ncchar('O');
        cdfobj.long_name = ncchar('OXYGEN (ML/L)            ');
        cdfobj.generic_name = ncchar('ox');
        cdfobj.FORTRAN_format = ncchar('f10.2');
        cdfobj.units = ncchar('ml/l');
        cdfobj.epic_code = nclong(60);
        cdfobj.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
        cdfobj.serial_number = ncchar(aa.chan8.serial);
        cdfobj.valid_range = ncfloat([0 5000]);
        cdfobj.cal_A = ncdouble(aa.chan8.cals(1));
        cdfobj.cal_B = ncdouble(aa.chan8.cals(2));
        cdfobj.cal_C = ncdouble(aa.chan8.cals(3));
        cdfobj.cal_D = ncdouble(aa.chan8.cals(4));
        cdfobj.salinity_setting = ncdouble(aa.chan8.salinity);
    elseif strcmp('%',aa.chan8.units) | strcmp('percent',aa.chan8.units),
        % % 62:OST:OXYGEN, %SAT             :ox:%: :
        cdf{'O_62'} = ncfloat('time', 'depth', 'lat', 'lon');
        cdf{'O_62'}.name = ncchar('OST');
        cdf{'O_62'}.long_name = ncchar('OXYGEN, %SAT          ');
        cdf{'O_62'}.generic_name = ncchar('ox');
        cdf{'O_62'}.FORTRAN_format = ncchar('f10.2');
        cdf{'O_62'}.units = ncchar('%');
        cdf{'O_62'}.epic_code = nclong(62);
        cdf{'O_62'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
        cdf{'O_62'}.serial_number = ncchar(aa.chan8.serial);
        cdf{'O_62'}.valid_range = ncfloat([0 1000]);
        cdf{'O_62'}.cal_A = ncdouble(aa.chan8.cals(1));
        cdf{'O_62'}.cal_B = ncdouble(aa.chan8.cals(2));
        cdf{'O_62'}.cal_C = ncdouble(aa.chan8.cals(3));
        cdf{'O_62'}.cal_D = ncdouble(aa.chan8.cals(4));
    else 
        disp('Warning: no recoginized o2 units, only raw data written to file')
    end

    % raw 02.  No EPIC key
    cdf{'o2raw'} = ncfloat('time', 'depth', 'lat', 'lon');
    cdf{'o2raw'}.name = ncchar('o2raw');
    cdf{'o2raw'}.long_name = ncchar('Raw Optode Output');
    cdf{'o2raw'}.generic_name = ncchar('o2raw');
    cdf{'o2raw'}.FORTRAN_format = ncchar('f10.2');
    cdf{'o2raw'}.units = ncchar('counts');
    cdf{'o2raw'}.epic_code = nclong(0);
    cdf{'o2raw'}.sensor_depth = ncdouble(metadata.water_depth-aa.instheight);
    cdf{'o2raw'}.serial_number = ncchar(aa.chan8.serial);
    cdf{'o2raw'}.valid_range = ncfloat([0 100]);
    cdf{'o2raw'}.salinity_setting = ncdouble(aa.chan8.salinity);
end

return

function add_vardesc(cdf)
% construct the EPIC attribute var_desc on the fly
theVars = var(cdf);
theDims = dim(cdf);
var_desc = '';
for ivar = 1:length(theVars),
    adim = 0;
    for idim = 1:length(theDims),
        if strcmp(ncnames(theVars{ivar}),ncnames(theDims{idim}))...
                | strcmp(ncnames(theVars{ivar}),'time2'),
            adim = 1;
        end
    end
    if ~adim,
        var_desc = sprintf('%s:%s',var_desc,theVars{ivar}.name(:));
    end
end
var_desc = var_desc(2:end);
cdf.VAR_DESC = ncchar(var_desc);
disp(sprintf('Writing variables %s',var_desc));
return

function add_fillvalues(cdf,fill)
theVars = var(cdf);
for i = 1:length(theVars),
    switch datatype(theVars{i}),
        case 'float'
            theVars{i}.FillValue_ = ncfloat(fill);
        case 'double'
            theVars{i}.FillValue_ = ncdouble(fill);
        case 'long'
            theVars{i}.FillValue_ = nclong(fill);
    end
end
return

% ---------------------------------------------
function add_minmaxvalues(cdf)
theVars = var(cdf);
for i = 1:length(theVars),
    if ~strcmp(ncnames(theVars{i}),'time') && ~strcmp(ncnames(theVars{i}),'time2'),
        data = theVars{i}(:);
        [row, col] = size(data);
        if col == 1,
            theVars{i}.minimum = ncfloat(min(data));
            theVars{i}.maximum = ncfloat(max(data));
        elseif col == 2,
            theVars{i}.minimum = ncfloat(min(min(data)));
            theVars{i}.maximum = ncfloat(max(max(data)));
        elseif col > 2, % then these are amp and cor
            for icol = 1:col,
                mins(icol) = min(data(:,icol));
                maxs(icol) = max(data(:,icol));
            end
            theVars{i}.minimum = ncfloat(mins);
            theVars{i}.maximum = ncfloat(maxs);
        end
    end
end
return

% ---------------------------------------------
function tjout = fix_time_nans(tjin,interval);
% tjin = time in Julian
% interval = interval between samples in sec
tjout = tjin;
dtjd = interval./(3600*24); % interval in julian day
badtj = find(isnan(tjin));
if length(badtj)  == 0, return; end
disp(sprintf('There are %d blank times to fill in',length(badtj)))
if badtj(1) == 1, % the first record is bad, this is a special case
    % find the first good time, then work backwards
    idx = 1;
    while badtj(idx) == idx,
        idx = idx+1;
    end
    idx = idx-1; % the last increment is not what we want, go back one
    % now idx = the last bad time in the beginning of the file
    for j = idx:-1:1
        tjout(j) = tjout(j+1)-dtjd;
    end
    % scan again for bad times
    badtj = find(isnan(tj));
    if length(bagtj == 0), return; end
end
for idx = 1:length(badtj),
    tjout(badtj(idx)) = tjout(badtj(idx)-1)+dtjd;    
end
return