function nc_loaddap_odv_surf
% Takes data from several surface instruments and puts it into an ASCII
% file that LOBOVIZ (or ODV) can read.  Currently takes the surface CTD,
% CO2, metsys, and chl data from netcdf files in DODS format and isus data
% from a text file.  
%
% NOTE: This program does not automatically find new deployments or find
% changes in NetCDF filenames.  Each change needs to be added to the
% load_data subroutine.
%
% Seth M. Bushinsky - 2007 - Modified 3/12/2008

urlbase = 'http://dods.mbari.org/cgi-bin/nph-nc/data/ssdsdata/deployments/m0/';
instrument_name = {'CTD0m','CO2QC','metsys','ocr3000'};%,'ecoflntu'};%'imCTD', 'ts';'ADCP-QC';'ecofl1';'ecofl2';'ecobb';'ocr3000';'gps';'ADCP';;},...

addpath('/usr/local/bin');
addpath('/ssdsdata/deployments/m0/ascii_all_dep/mfiles');

[cruise, data,  header, time] = nc_loaddap (urlbase, instrument_name);
%"-3" prevents the last few hours of data from being written because they often are NaNs as a
%result of interpolation.  
% 2007-04-19 Changes to "-13" because CO2 processing time is 12 hours off, so the
% last 1/2 day was being written as NaNs.  This lead to daily gaps in data.
%  Should change back to -3 once CO2 processing is fixed.
datalines = length(time)-13;
data = data(1:datalines,:);
cruise = cruise(1:datalines,:);

fidcfg_old = fopen('/ssdsdata/deployments/m0/ascii_all_dep/M0SURF.cfg', 'r');
if fidcfg_old ~= -1,
    A = fread(fidcfg_old, '*char');
    fclose(fidcfg_old);
    A = A';
    A = A(3:end);    
    %only append the new data to the text files.  
    prev_datalines = str2double(A);
    new_data = data(prev_datalines + 1 : datalines,:);
    new_cruise = cruise(prev_datalines + 1 : datalines,:);
    data = new_data;
    cruise = new_cruise;
end
    
if isempty(data),
    return,
end

%fidcfg = fopen('C:\Documents and Settings\sethb\My Documents\BOG\matlab\M0SURFtest.cfg', 'wt');
fidcfg = fopen('/ssdsdata/deployments/m0/ascii_all_dep/M0SURF.cfg', 'wt');
fprintf(fidcfg, '%s', ['//' int2str(datalines)]);
fclose(fidcfg);

if fidcfg_old ~= -1,
    fid = fopen('/ssdsdata/deployments/m0/ascii_all_dep/M0SURF.txt', 'r+');
    fprintf(fid, '%s', ['//' int2str(datalines)]);
else
    fid = fopen('/ssdsdata/deployments/m0/ascii_all_dep/M0SURF.txt', 'wt');
    fprintf(fid, '%s\n', ['//' int2str(datalines)]);
    for h = 1:length(header),
        if h~= length(header),
            fprintf(fid, '%s\t', char(header{h}));
        else
            fprintf(fid, '%s\n', char(header{h}));
        end
    end
end

fclose(fid);

%dlmwrite('C:\Documents and Settings\sethb\My Documents\BOG\matlab\M0SURFtest2.txt', data, 'delimiter','\t','-append', 'roffset', 2,'precision','%.4f');

disp('Calling tblwrite_smb() to write data to /ssdsdata/deployments/m0/ascii_all_dep/M0SURF.txt');
tblwrite_smb(data, header, cruise, '/ssdsdata/deployments/m0/ascii_all_dep/M0SURF.txt','\t',4,1)

end

%% Data Arrays

function [cruise, data,  header, time] = nc_loaddap (urlbase, instrument_name)
clear global time_hourly data_hourly
global time_hourly data_hourly

for i = 1:length(instrument_name),
    inst = char(instrument_name(i));

    %using the name of the instrument (inst) and the url base, load_data
    %retrieves the variable names and data contained in the netcdf files.
    %Data is obtained using loaddap and returned in a structure array
    %(inst_data).
    [variable_names, inst_data] = load_data (inst, urlbase);
    
    %loops through the multiple netcdf files for each instrument
    for f = 1:length(inst_data.(inst)),
        
        %this accomodates the differences between the quality controlled
        %netcdf file and the current netcdf file.  the current file has a
        %lot of entries that aren't used, so this removes them from the
        %data and variable name arrays.  
        if strcmp(char(instrument_name(i)), 'CO2QC') && f==2,
            for v = 1:length(variable_names.(inst){f}),
                attribute = char(variable_names.(inst){f}(v));
                
                if v>4 && v<49,
                    inst_data.(inst){f} = rmfield(inst_data.(inst){f}, (attribute));
                end
            end
            variable_names.(inst){f}(5:48)=[];
        end
        
        %loops through each field name for each netcdf file
        for v = 1:length(variable_names.(inst){f}),

            % This allows only the first data structure in ocr3000
            % (Eds_mean_410) to be read in.  Otherwise the script crashes
            % when the derived chl variable is introduced.  Also, no need
            % for all 40 or so variables to be read.
%             if strcmp(char(instrument_name(i)), 'ocr3000') && v>5,
%                continue, 
%             end   
            
            attribute = char(variable_names.(inst){f}(v));

            %checks that the attribute returns a structure, this indicates
            %that it is a multiple dimension variable and contains data
            attribute_struct = inst_data.(inst){f}.(attribute);
            if ~isstruct(attribute_struct),
                continue,
            end
            
            %puts the different fields within the attribute structure into
            %a cell array.  The first contains data, the second contains
            %time markers
            attribute_names = fieldnames(attribute_struct);
            
            %reads data into a cell/structure array for each attribute.
            %concatenates the data from each attribute across netcdf files.
            %Checks the size to determine if there are multiple data points
            %for each time, this indicates multiple depths.  
            if f==1,
                time_esecsx{i}.(attribute) = inst_data.(inst){f}.(attribute).(char(attribute_names(2)));
                temp = squeeze(inst_data.(inst){f}.(attribute).(char(attribute_names(1))));

                [m, n] = size(temp);

                if n == length(inst_data.(inst){f}.(attribute).(char(attribute_names(2)))),
                    temp = temp';
                end
                out_data{i}.(attribute) = temp;
            else
                time_esecsx{i}.(attribute) = [time_esecsx{i}.(attribute); inst_data.(inst){f}.(attribute).(char(attribute_names(2)))];
                temp = squeeze(inst_data.(inst){f}.(attribute).(char(attribute_names(1))));

                [m, n] = size(temp);

                if n == length(inst_data.(inst){f}.(attribute).(char(attribute_names(2)))),
                    temp = temp';
                end
                out_data{i}.(attribute) = [out_data{i}.(attribute); temp];
            end

        end
    end
    
        
    %interpolates data to an hourly grid
    disp(['Calling time_grid() for instrument ' inst]);
    time_grid (time_esecsx, out_data, i, inst);
    
    
    
end

% Add isus nitrate data to the gridded data
[isus_data isus_esecs isus_quality] = import_isus;
out_data{i+1}.('no3') = isus_data;
out_data{i+1}.('qual') = isus_quality;
time_esecsx{i+1}.('no3') = isus_esecs;
time_esecsx{i+1}.('qual') = isus_esecs;
time_grid (time_esecsx, out_data, i+1, 'isus');

        
time_edays = time_hourly/(3600*24);
epoch_start = datenum('1-jan-1970');
time = epoch_start + time_edays;

mdy = datestr(time, 'mm/dd/yyyy');
hhmm = datestr(time, 'HH:MM');

for row = 1:length(time),
    cruise{row,1} = 'M0';
    cruise{row,2} = num2str(row);
    cruise{row,3} = '*';
    cruise{row,4} = mdy(row,1:10);
    cruise{row,5} = hhmm(row,1:5);
end

lon(1:length(time),1) = inst_data.(inst){f}.(char(attribute_names(4)));
lat(1:length(time),1) = inst_data.(inst){f}.(char(attribute_names(3)));
%sets "Bot.Depth" to 0, purely for file structure benefit
depth(1:length(time),1) = 0;
days_2003 = time - (datenum('1-Jan-2003')) + 1;

data(:,1) = lon+360;
data(:,2) = lat;
data(:,3) = depth;
data(:,4) = days_2003;
data(:,5) = depth; %for now sensor depth is all 0, since this is the surface file... 

header = {'Cruise', 'Station', 'Type', 'mon/day/yr', 'hh:mm',...
    'Longitude [degrees_east]', 'Latitude [degrees_north]', 'Bot.Depth [m]',...
    'DaysSince2003[d]','SensorDepth[m]','QF', 'Temperature[°C]', 'QF', 'Salinity', 'QF', ...
    'AirPressure', 'QF', 'AirTemperature[°C]','QF',...
    'WindSpd[m/s]', 'QF', 'WindDir[°]', 'QF', ...
    'Nitrate[µM]', 'QF', 'Chlorophyll[µg/L]', 'QF',...
    'SeaPCO2[ppm]', 'QF', 'AirPCO2[ppm]', 'QF','OxygenSat[%]', 'QF', 'SolarRad[W/m2]', 'QF'};

data(:,6) = 0;

data(:,7) = data_hourly.CTD0m.Temperature;
data(:,8) = 0;
data(:,9) = data_hourly.CTD0m.Salinity;
data(:,10) = 0;
data(:,11) = data_hourly.metsys.AirPressure;
data(:,12) = 0;
miss_val = find(data(:,11)<-6000);
data(miss_val,11) = NaN;
data(:,13) = data_hourly.metsys.AirTemperature;
data(:,14) = 0;
data(:,15) = data_hourly.metsys.WindSpd_Windbird;
data(:,16) = 0;
data(:,17) = data_hourly.metsys.WindDir_Windbird;
data(:,18) = 0;
data(:,19) = data_hourly.isus.no3; %nitrate
data(:,20) = data_hourly.isus.qual;
data(:,21) = data_hourly.ocr3000.chl_490_Kd_0_10_mean;
data(:,22) = 0;
miss_val = find(data(:,21)<-1000);
data(miss_val,21) = NaN;
data(:,23) = data_hourly.CO2QC.sea_pCO2;
data(:,24) = 0;
miss_val = find(data(:,23)<0);
data(miss_val,23) = NaN;
data(:,25) = data_hourly.CO2QC.air_pCO2;
data(:,26) = 0;
miss_val = find(data(:,25)<0);
data(miss_val,25) = NaN;
data(:,27) = data_hourly.CO2QC.o2_sat;
data(:,28) = 0;
miss_val = find(data(:,27)<0);
data(miss_val,27) = NaN;
data(:,29) = data_hourly.ocr3000.Eds_mean_490; %Eds_mean_410
data(:,30) = 0;

end

%% Load Variables 
%% EDIT THIS ON EACH DEPLOYMENT
function [variable_names, instrument] = load_data(inst, urlbase)

switch char(inst);

    case {'CTD0m'}
        file.(inst){1} = '200405/m0_ctd0000_20040727_original.nc';
        file.(inst){2} = '200410/m0_ctd0000_20041006_original.nc';
        file.(inst){3} = '200504/m0_ctd0001_20050418_original.nc';
        file.(inst){4} = '200607/m0_ctd0001_20060731_original.nc';
        file.(inst){5} = '200701/m0_ctd0000_20070130_original.nc';
        file.(inst){6} = '200706/m0_ctd0001_20070621_original.nc';
    case {'imCTD'}
        file.(inst){1} = '200405/m0_imctd_20040526_original.nc';
        file.(inst){2} = '200410/m0_imctd_20041006_original.nc';
        file.(inst){3} = '200504/m0_imctd_20050418_original.nc';
        file.(inst){4} = '200607/m0_imctd_20060727_original.nc';
        file.(inst){5} = '200701/m0_imctd_20070130_original.nc';
        file.(inst){6} = '200706/m0_imctd_20070621_original.nc';
    case{'hourly'}
        file.(inst){1} = '200405/hourlyM0_20040525.nc';
        file.(inst){2} = '200410/hourlyM0_20041006.nc';
        file.(inst){3} = '200504/hourlyM0_20050418.nc';
        file.(inst){4} = '200607/hourlyM0_20060731.nc';
        file.(inst){5} = '200701/hourlyM0_20070130.nc';
        file.(inst){6} = '200706/hourlyM0_.nc';
    case{'metsys'}
        file.(inst){1} = '200405/m0_metsys_20040525_original.nc';
        file.(inst){2} = '200410/m0_metsys_20041006_original.nc';
        file.(inst){3} = '200504/m0_metsys_20050418_original.nc';
        file.(inst){4} = '200607/m0_metsys_20060801_original.nc';
        file.(inst){5} = '200607/m0_metsys_20061108_original.nc';
        file.(inst){6} = '200701/m0_metsys_20070130_original.nc';
        file.(inst){7} = '200706/m0_metsys_20070621_original.nc';
        file.(inst){8} = '200706/m0_metsys_20070716_original.nc';
    case{'CO2QC'}
        file.(inst){1} = 'qc/m0_pco2_qc.nc';
        file.(inst){2} = '200701/m0_pco2_20070130.nc';
        %mike mccann hasn't updated yet...
    case{'CO2'}
        file.(inst){1} = '200607/m0_pco2_20060731.nc';
        file.(inst){2} = '200701/m0_pco2_20070130.nc';
    case{'ADCP'}
        file.(inst){1} = '200405/m0_adcp1267_20040525.nc';
        file.(inst){2} = '200410/m0_adcp1267_20041006.nc';
        file.(inst){3} = '200504/m0_adcp1406_20050418.nc';
        file.(inst){4} = '200607/m0_adcp1267_20060731.nc';
        file.(inst){5} = '200701/m0_adcp1267_20070130.nc';
    case{'ecoflntu'}
        file.(inst){1} = '200607/m0_ecoflntu_20060731.nc';
        file.(inst){2} = '200607/m0_ecoflntu_20061108.nc';
        file.(inst){3} = '200701/m0_ecoflntu_20070130.nc';
        file.(inst){4} = '200706/m0_ecoflntu_20070621.nc';
    case{'ts'}
        file.(inst){1} = '200607/m0_ts_20060731.nc';
        file.(inst){2} = '200701/m0_ts_20070130.nc';
        file.(inst){3} = '200706/m0_ts_.nc';
    case{'ocr3000'}
        file.(inst){1} = '200405/m0_ocr3000_prev.nc';
        file.(inst){2} = '200607/m0_ocr3000_20060721.nc';
        file.(inst){3} = '200701/m0_ocr3000_20070130.nc';
        file.(inst){4} = '200706/m0_ocr3000_20070621.nc';
    case{'ecofl1'}
        file.(inst){1} = '200405/m0_ecofl_20040526.nc';
    case{'ecofl2'}
        file.(inst){1} = '200410/m0_ecofl_20041006.nc';
        file.(inst){2} = '200504/m0_ecofl_20050418.nc';
    case{'ecobb'}
        file.(inst){1} = '200504/m0_ecobb_20050418.nc';
        file.(inst){2} = '200504/m0_ecobb_20050728.nc';
    case{'gps'}
        file.(inst){1} = '200405/m0_gps1313_20040525.nc';
        file.(inst){2} = '200410/m0_gps1313_20041006.nc';
        file.(inst){3} = '200504/m0_gps1313_20050418.nc';
        file.(inst){4} = '200607/m0_gps1564_20060731.nc';
        file.(inst){5} = '200701/m0_gps1564_20070130.nc';
        file.(inst){6} = '200706/m0_gps1313_20070621.nc';
    case{'ADCP-QC'}
        file.(inst){1} = 'qc/200406/m0_adcp_20040603_workhorse.nc';
    case{'metsys-QC'}
        file.(inst){1} = 'qc/200406/m0_metsys_20040603.nc';
    case{'CTD-QC'}
        file.(inst){1} = 'qc/200406/m0_tstring_20040603.nc';
    case{'chl'}
        file.(inst){1} = '../../ocr3000cimt/netcdf/200701/chl_490_Kd_0_10_mean.nc';
        file.(inst){2} = '../../ocr3000cimt/netcdf/200706/chl_490_Kd_0_10_mean.nc';
end

for f = 1:length(file.(inst)),
    variable_names.(inst){f} = fieldnames(loaddap([urlbase char(file.(inst){f})]));
    instrument.(inst){f} = loaddap([urlbase char(file.(inst){f})]);
end
end

%% time gridding

%retrieves the names of different data types that are now stored in the
%out_data cell/structure array

function  time_grid (time_esecsx, out_data, i, inst)
data_names = fieldnames(out_data{i});

global time_hourly data_hourly

for d = 1:length(data_names),
    disp(['\nProcessing data_name ' char(data_names(d)) ':']); 
    time_esecs_start = time_esecsx{i}.(char(data_names(d)))(1);
    time_esecs_end = time_esecsx{i}.(char(data_names(d)))(end);

    %moves the start and end time to an even hour
    time_esecs_start = double((int32(time_esecs_start/3600)).*3600);
    time_esecs_end = double((int32(time_esecs_end/3600)+1).*3600);%the plus one establishes the end time after the last timestamp
    
    %sets the interval to which the data will be interpolated, in this
    %case 1 hour (60 min x 60 secs/min)
    interval = 3600;

    %sets the hourly grid to the start and end points of the first variable
    %it encounters.  Perhaps not the best plan, but a starting point.
    temp = time_esecs_start:interval:time_esecs_end;
    if isempty(time_hourly),
        time_hourly = temp;
    end
    
    %for some reason the time stamps aren't monotonic- possibly due to
    %overlap in deployments? this fixes that.
    %%y = diff([0;time_esecsx{i}.(char(data_names(d)))]);
    %%w = find(y>=0);
    %%mono_time = time_esecsx{i}.(char(data_names(d)))(w);
    %%mono_data = out_data{i}.(char(data_names(d)))(w);

    % Loop through all times, detect non-monotonic sections and save indices
    % into bindx (bad index) vector
    disp([ 'Finding non-monotonic sections in ' char(data_names(d))] )
    lastGoodTime = time_esecsx{i}.(char(data_names(d)))(1);
    bindx = [];
    for ii = 1:length(time_esecsx{i}.(char(data_names(d)))),
	if time_esecsx{i}.(char(data_names(d)))(ii) <= lastGoodTime,
		bindx = [bindx; ii];
	else
		lastGoodTime = time_esecsx{i}.(char(data_names(d)))(ii);
	end
    end
    disp([ 'Removing non-monotonic sections, length(bindx) = ' num2str(length(bindx))] )
    mono_time = time_esecsx{i}.(char(data_names(d)));
    mono_data = out_data{i}.(char(data_names(d)));
    mono_time(bindx) = [];
    mono_data(bindx) = [];
    
    if ~isempty(strfind(char(data_names(d)), 'Ed5_mean_410'))
	keyboard
    end
   
    % Iterate on removing non-monotonic time until non are left
    % This has the danger of taking forever.  Need to think about doing this better... (mpm 22 Sep 2008)
    neg = 1;
    while (0),
    %%while (isempty(neg)==0),
        q = diff([0;mono_time]);
        w = find(q>0);
        neg = find(q<=0);
	disp([ 'length(q) = ' num2str(length(q))  '  length(neg) = ' num2str(length(neg)) char(13)] )
        mono_time = mono_time(w);
        mono_data = mono_data(w);
    end
   
    disp(['Interpolating with interp1() data_name ' char(data_names(d)) ' (' num2str(d) ' of ' num2str(length(data_names)) ')']); 
    data_hourly.(inst).(char(data_names(d))) = interp1(mono_time, mono_data, time_hourly);

end
end

%% Import ISUS Nitrate data from a text file

function [isus_data epoch_secs isus_quality] = import_isus

C = importdata('/ssdsdata/isuscimt/data/M0.txt', ',');

yyyyddd = C.data(:,1);
hr = C.data(:,2);
isus_data = C.data(:,7);
isus_quality = C.data(:,13);

dec_year = yyyyddd_to_decyear(yyyyddd);

matlab_time = decimal_year2_matlab_time(dec_year);

matlab_hr = zeros(length(hr),1);
for k = 1:length(hr)
    matlab_hr(k) = datenum([0 0 0 hr(k) 0 0]);
end
   
matlab_time = matlab_time + matlab_hr;
epoch_secs = (matlab_time - datenum('1-jan-1970'))*3600*24;


end
